/* $Log: graph.c,v $
 * Revision 1.4  1995/03/20  14:54:55  martenl
 * *** empty log message ***
 *
 * Revision 1.3  1995/03/10  12:41:16  martenl
 * *** empty log message ***
 *
 * Revision 1.2  1995/02/23  09:45:18  martenl
 * getDrawBox:
 * 	Changed criteriom for when xmin,xmax etc is recomputed
 * 	to allow calling even if no mesh is present.
 *
 *
 */

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


#include "libsx.h"             /* gets us in the door with libsx          */
#include "globals.h"
#include "graph.h"
#include "grid.h"
#include "kraftwerk.h"
#include "util.h"
#include "naca.h"
#include "geodata.h"
#define min(x, y)                     (((x) < (y)) ? (x) : (y))
#define max(x, y)                     (((x) > (y)) ? (x) : (y))

#define PAPER_WIDTH 200       /* Paper size in mm */
#define PAPER_HEIGTH 300
#define PAPER_SIDE_MARGIN 0
#define PAPER_TOP_MARGIN 0
#define MM 2.8346

#define X11 0
#define PS 1

float   xmin_orig;
float   xmax_orig;
float   ymin_orig;
float   ymax_orig;

int device = X11;
float psPixScale;
float psPixXoff;
float psPixYoff;
float psPhysScale;
float psPhysXoff;
float psPhysYoff;

int cells[256];
int ncells = 0;
unsigned char red[255];
unsigned char green[255];
unsigned char blue[255];

FILE *psfile;
GC drawgc;
void psDrawFilledPolygon(XPoint *p, int n);

void
rainbow(double *hh, double *ss, double *vv, double *r, double *g, double *b);
void
dhsv2rgb(double h,double s, double v,double *r,double *g,double *b) ;

void myDrawLine(int x1, int y1, int x2, int y2)
{
  if(device == X11)
    DrawLine(x1, y1, x2, y2);
  else
    psDrawLine(x1, y1, x2, y2);
}
void myDrawCross(int x1, int y1, int h)
{
  if(device == X11){
    DrawLine(x1, y1-h, x1, y1+h);
    DrawLine(x1-h, y1, x1+h, y1);
  }
  else{
    psDrawLine(x1, y1-h, x1, y1+h);
    psDrawLine(x1-h, y1, x1+h, y1);
  }
}

void myDrawFilledBox(int x, int y, int w, int h)
{
  if(device == X11)
    DrawFilledBox(x, y, w, h);
  else
    psDrawFilledBox(x, y, w, h);
}
void myDrawFilledPolygon(XPoint *p, int n)
{
  if(device == X11)
    DrawFilledPolygon(p, n);
  else
    psDrawFilledPolygon(p, n);
}
void myDrawText(char *text, int x, int y)
{
  if(device == X11)
    DrawText(text, x, y);
  else
    psDrawText(text, x, y);
}
void drawText(char *text, float xp, float yp)
{
  int x, y;
  phys2pix(&x, &y, xp, yp);
  if(device == X11)
    DrawText(text, x, y);
  else
    psDrawText(text, x, y);
}

void line(float *x, float *y, int n)
{
  if(device == X11)
    xline(x, y, n);
  else
    psline(x, y, n);
}

void drawTitle(char *text)
{
  int x1, x2, y;
  int l;

  l = TextWidth(GetWidgetFont(drawWindow), text);
  phys2pix(&x1, &y, 0, ymax);
  phys2pix(&x2, &y, xmax, ymax);
  myDrawText(text, x1+(x2-x1-l)/2, y - 8);


}

/*
 * Coordinate transformation functions 
 *
 */
void pix2phys(int xx,int xy,float *px, float *py)
{
  float x, y;

  x= (float) (xx - Xxoff)/(float)Xscl;
  y= (float) (ymax+ymin) - (float) (xy - Xyoff)/(float)Xscl;

  *px = x;
  *py = y;
}


void phys2pix(int *x, int *y, float xx, float yy)
{
    *x = (int) ( xx * Xscl + Xxoff);
    *y = (int) ( (ymax+ymin-yy) * Xscl + Xyoff);
}


void pix2ps(int x, int y, float *xps, float *yps)
{
  *xps = (float) x * psPixScale + psPixXoff;
  *yps = PAPER_HEIGTH - PAPER_TOP_MARGIN - ((float) y * psPixScale + psPixYoff);
  *xps = *xps * MM;
  *yps = *yps * MM;

}


void phys2ps(float x, float y, float *xps, float *yps)
{
  int xpix, ypix;
  phys2pix(&xpix, &ypix, x, y);
  pix2ps(xpix, ypix, xps, yps);
}


void psPrint(char *printer, int where)
{
  char cmd[100];

  if(where == TO_PRINTER)
    psfile = fopen("plotBuf.ps", "w");
  else
    psfile = fopen(printer, "w");
  if(psfile != NULL){
    psInit();
    device = PS;
    doPlotting();
    psDone();
    fclose(psfile);

    if(where == TO_PRINTER){
      sprintf(cmd, "lpr -P%s plotBuf.ps", printer);
      system(cmd);
    }
    device = X11;
    doPlotting();
  }
  else
    errorHandler("Unable to open plotfile");
}


void psInit()
{

  psGetScale();
  fprintf(psfile, "%s\n", "%!ps-adobe-2.0");
  fprintf(psfile, "%s\n", "%%Creator: FEMLAB");
  fprintf(psfile, "%s\n", "%%Title: FEMLAB plot");
  fprintf(psfile, "%s\n", "%%Pages: 1");
  fprintf(psfile, "%s\n", "%%EndComments: 1");
  fprintf(psfile, "%s\n", "/mm { 2.8346 mul } def");

  fprintf(psfile, "%s\n", "/l { lineto } def");
  fprintf(psfile, "%s\n", "/m { moveto } def");
  fprintf(psfile, "%s\n", "/c { closepath } def");
  fprintf(psfile, "%s\n", "/n { newpath } def");
  fprintf(psfile, "%s\n", "/s { setgray } def");
  fprintf(psfile, "%s\n", "/r { setrgbcolor } def");
  fprintf(psfile, "%s\n", "/k { stroke } def");
  fprintf(psfile, "%s\n", "0.05 mm setlinewidth");

  fprintf(psfile, "%s\n", "/Helvetica-Bold findfont");
  fprintf(psfile, "%s\n", "30 scalefont");
  fprintf(psfile, "%s\n", "setfont");
  fprintf(psfile, "%s\n", "");

  fprintf(psfile, "%s\n", "0. setgray");
  fprintf(psfile, "%s\n", "76.5 mm 16.452 mm moveto");
  fprintf(psfile, "%s\n", "(F  E  M  L  A  B) show");
  fprintf(psfile, "%s\n", "5 mm 30 mm m");
  fprintf(psfile, "%s\n", "5 mm 295 mm l");
  fprintf(psfile, "%s\n", "200 mm 295 mm l");
  fprintf(psfile, "%s\n", "200 mm 30 mm l");
  fprintf(psfile, "%s\n", "5 mm 30 mm c");
  fprintf(psfile, "%s\n", "clip");
/*
  fprintf(psfile, "%s\n", "/Helvetica-Bold findfont");
  fprintf(psfile, "%s\n", "4.56 scalefont");
  fprintf(psfile, "%s\n", "setfont");
  fprintf(psfile, "%s\n", "77.6 mm 14 mm moveto");
  fprintf(psfile, "%s\n", "(T H E   C H O I C E  O F  A  N E W  G E N E R A T I O N) show");
*/
  fprintf(psfile, "%s\n", "newpath");

  fprintf(psfile, "%s\n", "/Courier findfont");
  fprintf(psfile, "%s\n", "10 scalefont");
  fprintf(psfile, "%s\n", "setfont");

}


void psDone()
{
  fprintf(psfile, "%s\n", "showpage");
}


void psGetScale()
{
  float scl;
  int width, height;
  float xoff, yoff;
  float ww;
  int nx, ny, id;


  if(state >= HAVE_SOLUTION){
    nx = nwinx;
    ny = nwiny;
    id = winid;
  }
  else{
    nx = ny = 1;
    id = 0;
  }

/* First fix the pixel -> paper scale */

  GetDrawAreaSize(&width, &height);
  width = width / nx;
  height = height / ny;

  scl  = (float) 1./max(width, height);
  ww = min((PAPER_WIDTH - 2*PAPER_SIDE_MARGIN)/nx, (PAPER_HEIGTH-2*PAPER_TOP_MARGIN)/ny);


  psPixScale  =  ww * scl;
  psPixXoff = 0*PAPER_SIDE_MARGIN;
  psPixYoff = 0*PAPER_TOP_MARGIN;

  psPixXoff += (PAPER_WIDTH/nx -  ( psPixScale*(width) + psPixXoff))/2  + (id+nx-1)%nx *PAPER_WIDTH/nx*scl;
  psPixYoff += (PAPER_HEIGTH/ny - ( psPixScale*(height) + psPixYoff))/2 + (int)(id/(nx+.5)) *PAPER_HEIGTH/ny*scl;

/* */
/* Now fix the physical -> paper scale */
/* */

  scl  = 1./max(xmax-xmin,ymax-ymin);
  xoff =  - xmin*scl;
  yoff =  - ymin*scl;

  psPhysScale =  ww * scl;
  psPhysXoff =  ww * xoff;
  psPhysYoff =  ww * yoff;
  psPhysXoff +=  (PAPER_WIDTH/nx -  ( psPhysScale*xmax +psPhysXoff))/2 + (id+nx-1)%nx *PAPER_WIDTH/nx;
  psPhysYoff +=  (PAPER_HEIGTH/ny - ( psPhysScale*ymax +psPhysYoff))/2 + (int)(id/(nx+.5)) *PAPER_HEIGTH/ny;

}

/*
 * Color stuff
 * 
 *
 */

mySetColor(int color)
{
    if(have_colors == TRUE){
      SetColor(color);
      if(color == WHITE){
	psred = 0.0;
	psgreen = 0.0;
	psblue = 0.0;
      }
      else if(color == BLACK){
	psred = 1.0;
	psgreen = 1.0;
	psblue = 1.0;
      }
      else if(color == RED){
	psred = 1.0;
	psgreen = 0.0;
	psblue = 0.0;
      }
      else if(color == GREEN){
	psred = 0.0;
	psgreen = 1.0;
	psblue = 0.0;
      }
      else if(color == BLUE){
	psred = 0.0;
	psgreen = 0.0;
	psblue = 1.0;
      }
      else if(color == YELLOW){
	psred = 0.0;
	psgreen = 1.0;
	psblue = 1.0;
      }
      else if(color == GREY){
	psred = .8;
	psgreen = .8;
	psblue = .8;
      }
    }

}
void mySetMyColor(float value)
{
  if(have_colors == TRUE){
    psred   = (float)red[(int) ((float)(value*ncells))]/255.;
    psgreen = (float)green[(int) ((float)(value*ncells))]/255.;
    psblue  = (float)blue[(int) ((float)(value*ncells))]/255.;
    SetColor(cells[(int) ((float)ncells * value)]);
  }
}





/*
 *   Postscript counterparts of plotting functions
 *
 */

void psDrawFilledBox(int x, int y, int w, int h)
{
  float xx, yy;

  fprintf(psfile,"%s","n\n");
  fprintf(psfile,"%s","0. s\n");
  pix2ps(x, y, &xx, &yy);
  fprintf(psfile,"%f %f m\n", xx, yy);
  pix2ps(x + h, y, &xx, &yy);
  fprintf(psfile,"%f %f l\n", xx, yy);
  pix2ps(x + h, y + h, &xx, &yy);
  fprintf(psfile,"%f %f l\n", xx, yy);
  pix2ps(x, y + h, &xx, &yy);
  fprintf(psfile,"%f %f l\n", xx, yy);
  pix2ps(x, y, &xx, &yy);
  fprintf(psfile,"%f %f l\n", xx, yy);

  fprintf(psfile,"%s","c\n");
  fprintf(psfile,"%s",".0 s fill\n");

}
void psDrawFilledPolygon(XPoint *p, int n)
{
  float x, y;
  int i;
  fprintf(psfile,"%s","n\n");
  if(color_ps == TRUE)
    fprintf(psfile,"%f %f %f r\n", psred, psgreen, psblue);
  else
    fprintf(psfile,"%s","0. s\n"); 

  pix2ps(p[0].x, p[0].y, &x, &y);
  fprintf(psfile,"%f %f m\n", x, y);
  for(i=1; i < n; i++){
    pix2ps(p[i].x, p[i].y, &x, &y);
    fprintf(psfile,"%f %f l\n", x, y);
  }
  fprintf(psfile,"c\n");
  fprintf(psfile,"fill\n");

}

void psDrawText(char *text, int x, int y)
{
  float xx, yy;

  fprintf(psfile,"%s","0. s\n");

  pix2ps(x, y, &xx, &yy);
  
  fprintf(psfile, "%f %f m\n", xx, yy);
  fprintf(psfile, "(%s) show \n", text);

}

void psDrawLine(int x1, int y1, int x2, int y2)
{
  float xx, yy;

  fprintf(psfile,"%s","n\n");
  fprintf(psfile,"%s","0. s\n");

  pix2ps(x1, y1, &xx, &yy);
  fprintf(psfile,"%.1f %.1f m\n", xx, yy);

  pix2ps(x2, y2, &xx, &yy);
  fprintf(psfile,"%.1f %.1f l\n", xx, yy);

  fprintf(psfile,"%s","k\n");
}
void psline(float *x, float *y,int np)
{
  float xx, yy;
  int i;

  fprintf(psfile,"%s","n\n");
  if(color_ps == TRUE)
    fprintf(psfile,"%f %f %f r\n", psred, psgreen, psblue);
  else
    fprintf(psfile,"%s","0. s\n"); 

  phys2ps(*x++, *y++, &xx, &yy);
  fprintf(psfile,"%.1f %.1f m\n", xx, yy);

  for(i = 1; i < np; i++){
    phys2ps(*x++, *y++, &xx, &yy);
    fprintf(psfile,"%.1f %.1f l\n", xx, yy);
  }
  fprintf(psfile,"%s","k\n");
}





void getDrawBox()
{
  int i, j, n, nb;
  PTS_NODE *p;

  if(zoomActive == FALSE && state <= HAVE_BOUNDARY){
    xmin = xmin_orig;
    xmax = xmax_orig;
    ymin = ymin_orig;
    ymax = ymax_orig;
  }

  if(zoomActive == FALSE && state >= HAVE_MESH){
    xmin = 10000;
    ymin = 10000;
    xmax = -10000;
    ymax = -10000;
    FOR_ALL_POINTS(p){
      xmin = min(xmin, p->physicalX);
      ymin = min(ymin, p->physicalY);
      xmax = max(xmax, p->physicalX);
      ymax = max(ymax, p->physicalY);
    }      
    if(zoomActive == FALSE && state >= HAVE_SOLUTION){
      xmin = 10000;
      ymin = 10000;
      xmax = -10000;
      ymax = -10000;
      for(i=0; i < mySol[currentSolution].nnode; i++){
	xmin = min(xmin, mySol[currentSolution].coord[i][0]);
	ymin = min(ymin, mySol[currentSolution].coord[i][1]);
	xmax = max(xmax, mySol[currentSolution].coord[i][0]);
	ymax = max(ymax, mySol[currentSolution].coord[i][1]);
      }
    }
  }
  
}


void xline(float *x,float *y, int np)
{
  int i;
  XPoint *pp, *pts;

  pts = (XPoint *) malloc( np * sizeof(XPoint) );
  pp = pts;

  for(i = 0; i < np; i++)  {
    pp[i].x = (int) ( *x++ * Xscl + Xxoff);
    pp[i].y = (int) ( (ymax+ymin-*y++) * Xscl + Xyoff);
  }
  DrawPolyline(pp, np);
  free(pts);
} 
void fill(float *x,float *y, int np)
{
  int i;
  XPoint *pp, *pts;

  pts = (XPoint *) malloc( np * sizeof(XPoint) );
  pp = pts;

  for(i = 0; i < np; i++)  {
    pp[i].x = (int) ( *x++ * Xscl + Xxoff);
    pp[i].y = (int) ( (ymax+ymin-*y++) * Xscl + Xyoff);
  }
  DrawFilledPolygon(pp, np);
  free(pts);
} 
void arrow(float xr1, float yr1, float xr2, float yr2)
{
  float x[5], y[5];
  double dx, dy, d1, d;
    
  x[0]=xr2;
  x[1]=xr1;
  y[0]=yr2;
  y[1]=yr1;
  line(x,y,2);

  dx=xr2-xr1;                                                              
  dy=yr2-yr1;
  d1=sqrt(dx*dx+dy*dy);
  d = .1 * d1;
  if(d1 != 0.){
    dx=dx/d1;
    dy=dy/d1;
    x[1]=xr2+(-dy-dx)*d;
    y[1]=yr2+(dx-dy)*d;
    x[2]=xr2;
    y[2]=yr2;
    x[3]=xr2+(dy-dx)*d;
    y[3]=yr2+(-dx-dy)*d;
    x[4]=xr2;
    y[4]=yr2;
    line(x,y,5);
  }

}

void setClipRect(int x1, int y1, int x2, int y2)
{
  XRectangle box[1];
  Display        *dpy;
  int fore_g,back_g;      /* Fore and back ground pixels */
  GC  drawgc;
 
  dpy = (Display *) XtDisplay(drawWindow);

  back_g = WhitePixel(XtDisplay(drawWindow),DefaultScreen(XtDisplay(drawWindow)));
  fore_g = BlackPixel(XtDisplay(drawWindow),DefaultScreen(XtDisplay(drawWindow)));
 
  /* Create drawing GC */
  drawgc = XCreateGC(XtDisplay(drawWindow), DefaultRootWindow(XtDisplay(drawWindow)), 0, 0);
 

  box[0].x = x1;
  box[0].y = y1;
  box[0].width = x2 - x1;
  box[0].height = y2 - y1;

  printf("%i %i %i %i\n", x1, y1, x2, y2);
  XSetClipRectangles(dpy, drawgc, 0, 0, box, 1, Unsorted);


}


void getDrawingScale()
{
  float xoff, yoff;
  float iscl;	/* inverse scale */
  int width, height;
  int ww;
  int nx, ny, id;

  if(state >= HAVE_MESH)
    getDrawBox();

  if(state >= HAVE_SOLUTION){
    nx = nwinx;
    ny = nwiny;
    id = winid;
  }
  else{
    nx = ny = 1;
    id = 0;
  }
  GetDrawAreaSize(&width, &height);

  width = width / nx;
  height = height / ny;

  iscl  = max(xmax-xmin,ymax-ymin);
  if (iscl != 0.0) {
   if(state > HAVE_NOTHING) {
/*  These are recalculated ??
    xoff = .5-(xmin+.5*(xmax-xmin))/iscl;
    yoff = .5-(ymin+.5*(ymax-ymin))/iscl;
*/
      xoff =  - xmin*.85/iscl;
      yoff =  - ymin*.85/iscl;
      ww = min(width, height);
      Xscl  = (int) ((float)ww * 0.85 / iscl);
      Xxoff = (int) ww * xoff ;
      Xyoff = (int) ww * yoff;
      Xxoff += (int) (width -  ( Xscl*xmax +Xxoff))/2 + (id+nx-1)%nx *width;
      Xyoff += (int) (height - ( Xscl*ymax +Xyoff))/2 + (int)(id/(nx+.5)) *height;
   }

  SetDrawArea(drawWindow);
  if(device == PS)
    psGetScale();
  }
}

void getBoundaryScale(Widget w, void *data)
{
  int ww=0, hh=0;

  TagList tags[] = {
    {TAG_WINDOW_LABEL,  "Define the size of your domain ", NULL    , TAG_INIT}, 
    {TAG_LABEL,  "Before the boundary can be defined",     NULL,  TAG_NOINIT},
    {TAG_LABEL,  "You must specify the size of your domain",     NULL,  TAG_NOINIT},
    {TAG_FLOAT, "min X:   ", &xmin, TAG_INIT},
    {TAG_FLOAT, "max X:   ", &xmax, TAG_INIT},
    {TAG_FLOAT, "min Y:   ", &ymin, TAG_INIT}, 
    {TAG_FLOAT, "max Y:   ", &ymax, TAG_INIT}, 
    {TAG_DONE,   NULL      , NULL,  TAG_NOINIT}
  };
  xmax = 1.;
  ymax = 1.;
  xmin = 0.;
  ymin = 0.;

  if(GetValues(tags))
    printf("Cancelled\n");

  getDrawingScale();
  if(state == HAVE_NOTHING)
    state = HAVE_SCALE;

  xmin_orig = xmin;
  xmax_orig = xmax;
  ymin_orig = ymin;
  ymax_orig = ymax;
  redisplay(w, ww, hh, data);
  naca(w, data);
}



void freeMyColorMap(int n)
{
  int i;

  for(i = 0; i < n ; i++)
    FreePrivateColor(cells[i]);

}

int setUpMycolorMap(int ncolor)
{
  double h, s, v, r, g, b;
  int i, n;

  n = ncolor;
  freeMyColorMap(ncells);

/* Get all available colors */
  i = 0;
  while( (cells[i++] = GetPrivateColor()) != -1 );
  n = i - 1;

  ncells = n;

  h=1.;
  s=1.;
  v=1.;
  for(i = 0; i < n ; i++){
    h = 1.0 - 1.0* (float)i /( n );
/*    h = .5 + .5*h;*/
    rainbow(&h,&s,&v,&r,&g,&b);
    red[i]   = 255*r;
    green[i] = 255*g;
    blue[i]  = 255*b;
  }
  for(i = 0; i < n ; i++)
    SetPrivateColor(cells[i], red[i], green[i], blue[i]);

  return n;

}



/*   rainbow(h, s, v, r, g, b)
     double h, s, v, *r, *g, *b;

 This routine computes colors suitable for use in color level plots.
 Typically s=v=1 and h varies from 0 (red) to 1 (blue) in
 equally spaced steps.  (h=.5 gives green; 1<h<1.5 gives magenta.)
 To convert for frame buffer, use   R = floor(255.999*pow(*r,1/gamma))  etc.
 To get tables calibrated for other devices or to report complaints,
 contact ehg@research.att.com.
*/

/*
 * The author of this software is Eric Grosse.  Copyright (c) 1986,1991 by AT&T.
 * Permission to use, copy, modify, and distribute this software for any
 * purpose without fee is hereby granted, provided that this entire notice
 * is included in all copies of any software which is or includes a copy
 * or modification of this software and in all copies of the supporting
 * documentation for such software.
 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
 * WARRANTY.  IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY
 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
 * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
 */

double	huettab[] = {
 0.0000, 0.0062, 0.0130, 0.0202, 0.0280, 0.0365, 0.0457, 0.0559, 0.0671, 0.0796,
 0.0936, 0.1095, 0.1275, 0.1482, 0.1806, 0.2113, 0.2393, 0.2652, 0.2892, 0.3119,
 0.3333, 0.3556, 0.3815, 0.4129, 0.4526, 0.5060, 0.5296, 0.5501, 0.5679, 0.5834,
 0.5970, 0.6088, 0.6191, 0.6281, 0.6361, 0.6430, 0.6490, 0.6544, 0.6590, 0.6631,
 0.6667, 0.6713, 0.6763, 0.6815, 0.6873, 0.6937, 0.7009, 0.7092, 0.7190, 0.7308,
 0.7452, 0.7631, 0.7856, 0.8142, 0.8621, 0.9029, 0.9344, 0.9580, 0.9755, 0.9889,
 1.0000
};
/* computed from the FMC-1 color difference formula */
/* Barco monitor, max(r,g,b)=1, n=61 magenta,  2 Jan 1986 */

void rainbow(hh, ss, vv, r, g, b)
double	*hh, *ss, *vv, *r, *g, *b;
{
	int	i;
	double h,s,v;
	double	modf();
	double trash;
	h = *hh;
	s = *ss;
	v = *vv;

	h = 60 * modf(h / 1.5, &trash);
	i = floor(h);
	h = huettab[i] + (huettab[i+1] - huettab[i]) * (h - i);
	dhsv2rgb(h, s, v, r, g, b);
}

void
dhsv2rgb(h, s, v, r, g, b)    /*...hexcone model...*/
double	h, s, v, *r, *g, *b;    /* all variables in range [0,1[ */
/* here, h=.667 gives blue, h=0 or 1 gives red. */
{  /* see Alvy Ray Smith, Color Gamut Transform Pairs, SIGGRAPH '78 */
	int	i;
	double	f, m, n, k;
	double	modf();
	double trash;
	h = 6 * modf(h, &trash);
	i = floor(h);
	f = h - i;
	m = (1 - s);
	n = (1 - s * f);
	k = (1 - (s * (1 - f)));
	switch (i) {
	case 0: 
		*r = 1; 
		*g = k; 
		*b = m; 
		break;
	case 1: 
		*r = n; 
		*g = 1; 
		*b = m; 
		break;
	case 2: 
		*r = m; 
		*g = 1; 
		*b = k; 
		break;
	case 3: 
		*r = m; 
		*g = n; 
		*b = 1; 
		break;
	case 4: 
		*r = k; 
		*g = m; 
		*b = 1; 
		break;
	case 5: 
		*r = 1; 
		*g = m; 
		*b = n; 
		break;
	default: 
		fprintf(stderr, "bad i: %f %d", h, i); 
		exit(1);
	}
	f = *r;
	if ( f < *g ) 
		f = *g;
	if ( f < *b ) 
		f = *b;
	f = v / f;
	*r *= f;
	*g *= f;
	*b *= f;
}


syntax highlighted by Code2HTML, v. 0.9.1