/* psgraph.c      functions for plot data in postscript lang.
				1999  8 31  Tamagawa	*/
				/* 2000 6 7 mod.  axis */

#include "psgraph.h"


/*    FUNCTIONS  FOR PLOT IN POSTSCRIPT     */



		
FILE *  psinit( char file[])
{
	FILE *fp;
	fp = fopen(file, "w");

	psgraph.space = 1.0;		/* 1cm space */
	psgraph.papersize_x = 21.0;   		/* A4 */
	psgraph.papersize_y = 30.0; 
	psgraph.logx = psgraph.logy = 0;	/* default linear */



	/* set default */
	strcpy(psgraph.fontname, "Helvetica");
	psgraph.fontsize = 18;
	psgraph.factor_xa = 1 / 2.54 * 72;
	psgraph.factor_xb = 0;
	psgraph.factor_ya = 1 / 2.54 * 72;
	psgraph.factor_yb = 0;

	/* fprintf(fp, "%%!PS-Adobe-2.0\n"); */
	fprintf(fp, "/\%s findfont\n", psgraph.fontname);
	fprintf(fp, "%.2f scalefont\n", psgraph.fontsize);
	fprintf(fp, "setfont\n");

	return fp;

}

void psclose(FILE *fp)
{
	fprintf(fp, "showpage\n");
	fclose(fp);
}

void pssetlogx( int flag )
{
	if(flag ==1 || flag == 0 || flag == -1 ) {
		psgraph.logx = flag;
	} else {
		fprintf(stderr,"Wrong arg (pssetlogx)\n");
	}
}
void pssetlogy( int flag )
{
	if(flag ==1 || flag == 0 || flag == -1 ) {
		psgraph.logy = flag;
	} else {
		fprintf(stderr,"Wrong arg (pssetlogx)\n");
	}
}



void pssetwindow(float a_xmin, float a_xmax, float a_ymin, float a_ymax,
				float xmin, float xmax, float ymin, float ymax)
{
	/* set frame to draw   (a_xmin, a_ymin) - (a_xmax, a_xmin) ==>
							(xmin.ymin) - (xmax, ymax) */
	/*   xmin,xmax,ymin,ymax  in cm  */


	if( a_xmin <  psgraph.space  ) a_xmin = psgraph.space;
	if( a_xmax >  psgraph.papersize_x -  psgraph.space )
		a_xmax = psgraph.papersize_x - psgraph.space;
	if( a_ymin <  psgraph.space  ) a_ymin = psgraph.space;
	if( a_ymax >  psgraph.papersize_y -  psgraph.space )
		a_ymax = psgraph.papersize_y - psgraph.space;

	xmin = pslogch(xmin, psgraph.logx );
	xmax = pslogch(xmax, psgraph.logx );
	ymin = pslogch(ymin, psgraph.logy );
	ymax = pslogch(ymax, psgraph.logy );

	psgraph.factor_xa = (a_xmax - a_xmin)/(xmax - xmin) / 2.54 * 72;
	psgraph.factor_xb = a_xmin / 2.54 * 72 - xmin * psgraph.factor_xa;
	psgraph.factor_ya = (a_ymax - a_ymin)/(ymax - ymin) / 2.54 * 72;
	psgraph.factor_yb = a_ymin / 2.54 * 72 - ymin * psgraph.factor_ya;

}

float pslogch( float x, int flag )
{
	if( flag == 1 ) {
		return log10(x);
	} else if (flag == -1 ) {
		return log10(-x);
	} else {
		return x;
	}
}

void psgetpos( float x, float y, float *xcm, float *ycm)
{
	*xcm = (psgraph.factor_xa * pslogch(x, psgraph.logx)
					+ psgraph.factor_xb)/CM_IN_PT;
	*ycm = (psgraph.factor_ya * pslogch(y, psgraph.logy)
					+ psgraph.factor_yb)/CM_IN_PT;
}
	


void psmove( FILE *fp, float x, float y)
{
	float dx , dy;

	x = pslogch( x, psgraph.logx );
	y = pslogch( y, psgraph.logy );

	dx = psgraph.factor_xa * x + psgraph.factor_xb;
	dy = psgraph.factor_ya * y + psgraph.factor_yb;
	
	fprintf(fp, " %.2f %.2f moveto\n", dx, dy);
}

void psmove_cm( FILE *fp, float x, float y)
{
	float dx , dy;


	dx = CM_IN_PT * x;
	dy = CM_IN_PT * y;
	
	fprintf(fp, " %.2f %.2f moveto\n", dx, dy);
}

void psdraw( FILE *fp, float x, float y)
{
	float dx , dy;

	x = pslogch( x, psgraph.logx );
	y = pslogch( y, psgraph.logy );
	dx = psgraph.factor_xa * x + psgraph.factor_xb;
	dy = psgraph.factor_ya * y + psgraph.factor_yb;


	fprintf(fp, " %.2f %.2f lineto\n", dx, dy);
}

void psdraw_cm( FILE *fp, float x, float y)
{
	float dx , dy;

	dx = CM_IN_PT * x;
	dy = CM_IN_PT * y;

	fprintf(fp, " %.2f %.2f lineto\n", dx, dy);
}

void pscircle( float x, float y, float rad, FILE *fp )
{

	x = pslogch( x, psgraph.logx );
	y = pslogch( y, psgraph.logy );

	fprintf( fp, " %.2f %.2f  moveto  %.2f %.2f %.2f 0 360 arc fill\n",
		psgraph.factor_xa * x + psgraph.factor_xb,
		psgraph.factor_ya * y + psgraph.factor_yb,
		psgraph.factor_xa * x + psgraph.factor_xb,
		psgraph.factor_ya * y + psgraph.factor_yb,
		CM_IN_PT * rad );
}


void psbox(float xmin, float xmax, float ymin, float ymax, FILE * fp )
{
/* draw box 	*/

	fprintf( fp, "newpath\n");
	psmove( fp, xmin, ymin );
	psdraw( fp, xmin, ymax );
	psdraw( fp, xmax, ymax );
	psdraw( fp, xmax, ymin );
	psdraw( fp, xmin, ymin );
	fprintf( fp, "closepath\n");
	fprintf( fp, "stroke\n");
}


void psline( float *x, float *y, int num, FILE *fp)
{
	int i;

	fprintf( fp, "newpath\n");
	psmove( fp, *x, *y );
	for(i=0; i < num; i++) {
		psdraw( fp, *(x+i), *(y+i) );
	}
	fprintf( fp, "stroke\n");
}




void pslinewidth( int width, FILE *fp)
{

	fprintf( fp, "%d setlinewidth\n", width);
}

void pslinegray( int gray, FILE *fp)
{
	fprintf( fp, "%d setgray\n", gray);
}
	
void pssetfont( char fontname[], float scale, FILE * fp )
{
		
	strcpy(psgraph.fontname, fontname);
	psgraph.fontsize = scale;
	fprintf(fp, "/%s findfont\n", psgraph.fontname);
	fprintf(fp, "%.0f scalefont\n", psgraph.fontsize);
	fprintf(fp, "setfont\n");

}

void psnumberi( float x, float y, int number, char format[],
					int pos, int angle,  FILE *fp)
{
	char str[MAXSTRING];
	sprintf( str, format, number );
	psstrings( x, y, str, pos, angle, fp);
}
void psnumberf( float x, float y, float number, char format[], int pos,
		int angle, FILE *fp)
{
	char str[MAXSTRING];
	sprintf( str, format, number );
	psstrings( x, y, str, pos, angle, fp);
}

void psnumberi_cm( float x, float y, int number, char format[],
					int pos, int angle,  FILE *fp)
{
	char str[MAXSTRING];
	sprintf( str, format, number );
	psstrings_cm( x, y, str, pos, angle, fp);
}
void psnumberf_cm( float x, float y, float number, char format[], int pos,
		int angle, FILE *fp)
{
	char str[MAXSTRING];
	sprintf( str, format, number );
	psstrings_cm( x, y, str, pos, angle, fp);
}

void psstrings_cm( float x, float y, char strings[], int pos, int angle, FILE *fp)
{
/*     pos :         1      2      3
                     4      5      6
		     7      8      9		*/
/*  (x, y)  is at    "pos" of the string (cm !)	*/
/*   angle     degree.				*/



	fprintf( fp, "%.2f %.2f moveto\n", CM_IN_PT * x, CM_IN_PT * y);

	fprintf( fp, "gsave\n");

	fprintf( fp, "currentpoint  translate %d rotate\n", angle);

	if( pos == 1 || pos == 2 || pos == 3 )
		fprintf( fp, "0 %.2f rmoveto ", -psgraph.fontsize);
	if( pos == 4 || pos == 5 || pos == 6 )
		fprintf( fp, "0 %.2f rmoveto ", -0.5 * psgraph.fontsize);

	if( pos == 1 || pos == 4 || pos == 7 )
	     fprintf( fp, "(%s) show\n", strings);
	if( pos == 2 || pos == 5 || pos == 8 )
	     fprintf( fp, "(%s) dup stringwidth pop -0.5 mul 0 rmoveto show\n",
								strings);
	if( pos == 3 || pos == 6 || pos == 9 )
		fprintf( fp, "(%s) dup stringwidth pop -1 mul 0 rmoveto show\n",
								strings);


	fprintf( fp, "grestore\n");
}

void psstrings( float x, float y, char strings[], int pos, int angle, FILE *fp)
{
/*     pos :         1      2      3
                     4      5      6
		     7      8      9		*/
/*  (x, y)  is at    "pos" of the string	*/
/*   angle     degree.				*/


	psmove( fp, x, y);

	fprintf( fp, "gsave\n");

	fprintf( fp, "currentpoint  translate %d rotate\n", angle);

	if( pos == 1 || pos == 2 || pos == 3 )
		fprintf( fp, "0 %.2f rmoveto ", -psgraph.fontsize);
	if( pos == 4 || pos == 5 || pos == 6 )
		fprintf( fp, "0 %.2f rmoveto ", -0.5 * psgraph.fontsize);

	if( pos == 1 || pos == 4 || pos == 7 )
	     fprintf( fp, "(%s) show\n", strings);
	if( pos == 2 || pos == 5 || pos == 8 )
	     fprintf( fp, "(%s) dup stringwidth pop -0.5 mul 0 rmoveto show\n",
								strings);
	if( pos == 3 || pos == 6 || pos == 9 )
		fprintf( fp, "(%s) dup stringwidth pop -1 mul 0 rmoveto show\n",
								strings);


	fprintf( fp, "grestore\n");
}


void psxaxis( float xmin, float xmax, float dx, float y,
	float label_min, float label_max, float d_label, char form[], 
	char which[], FILE *fp)
{
/*   draw X - axis  from (xmin,y) to (amax, y)
	label number is written from label_min to label_max
	every d_label with format form[]			*/
/* which   =   b or t  or n    t means  top  axis, n means no label */
	float t, y1, y2;
	float ax_length = 5/CM_IN_PT ;		/*    5 pt	*/

	float xcm, ycm;
	int i;
	int labelpos, ls;
	float tref;


	fprintf( fp, "newpath\n");
	psgetpos(xmin, y, &xcm, &ycm );
	psmove_cm( fp, xcm, ycm );
	psgetpos(xmax, y, &xcm, &ycm );
	psdraw_cm( fp, xcm, ycm );

	if( strcmp( which, "t" ) == 0 || strcmp( which, "T") == 0) {
		labelpos = 8; ls = 1;
	} else {
		labelpos = 2; ls = -1;
	}

	if( psgraph.logx != 1 && psgraph.logx != -1 ) {
		if( dx >0 ) {
		for( t = xmin; t <= xmax; t += dx ) {
		psgetpos(t, y, &xcm, &ycm );
		psmove_cm( fp, xcm, ycm+ax_length );
		psdraw_cm( fp, xcm, ycm-ax_length );
		} 
		}
		fprintf( fp, "closepath stroke\n" );
		if(!( strcmp( which, "n" ) == 0  
				|| strcmp( which, "N") == 0 )) {
			if( d_label > 0 ) {
			for( t =label_min; t <= label_max; t += d_label) {
			psgetpos(t, y, &xcm, &ycm );
			psnumberf_cm( xcm, ycm + ls * 1.1 * ax_length,
						t, form, labelpos, 0, fp );
			}
			}
		}
	} else {   /* log scale */
		if( dx > 0 ) {
		for( tref = xmin, t = xmin; t <= xmax; t += dx ) {
		psgetpos(t, y, &xcm, &ycm );
		psmove_cm( fp, xcm, ycm+ax_length );
		psdraw_cm( fp, xcm, ycm-ax_length );
		if( t > tref*10 ) {
			dx = 10*dx;
			tref = t;
		}

		}
		}
		fprintf( fp, "closepath stroke\n" );
		if(!( strcmp( which, "n" ) == 0  
			|| strcmp( which, "N") == 0 )) {
		if( d_label > 0 ) {
		for( t =label_min; t <= label_max; t *= d_label) {
			psgetpos(t, y, &xcm, &ycm ); 
			psnumberf_cm( xcm, ycm + ls * 1.1 * ax_length,
						t, form, labelpos, 0, fp );
			}
		}
		}
	}

	
}


void psyaxis( float ymin, float ymax, float dy, float x,
	float label_min, float label_max, float d_label, char form[],
	char which[], FILE *fp)
{
/*   draw Y - axis  from (x, ymin) to (x, ymax)
	label number is written from label_min to label_max
	every d_label with format form[]			*/
/* which   =   l or r  or n    r means  Right axis, n means no label */

	float t, x1, x2;
	float ax_length = 5/CM_IN_PT ;          /*    5 pt      */

	float xcm, ycm;
	int i;
	int labelpos, ls;
	float tref;


	fprintf( fp, "newpath\n");
	psgetpos( x, ymin, &xcm, &ycm );
	psmove_cm( fp, xcm, ycm );
	psgetpos( x, ymax, &xcm, &ycm );
	psdraw_cm( fp, xcm, ycm );

	if( strcmp( which, "r" ) == 0  || strcmp( which, "R") == 0 ) {
		labelpos = 4; ls = 1;
	} else {
		labelpos = 6; ls = -1;
	}
		
	if( psgraph.logy != 1 && psgraph.logy != -1 ) {
		if( dy > 0 ) {
		for( t = ymin; t <= ymax; t += dy ) {
		psgetpos( x, t, &xcm, &ycm );
		psmove_cm( fp, xcm+ax_length, ycm );
		psdraw_cm( fp, xcm-ax_length, ycm );
		} 
		}
		fprintf( fp, "closepath stroke\n" );
		if(!( strcmp( which, "n" )==0 
			|| strcmp( which, "N")==0 )) {
			if( d_label > 0 ) {
			for( t =label_min; t <= label_max; t += d_label) {
			psgetpos( x, t, &xcm, &ycm );
			psnumberf_cm( xcm + ls * 1.3 * ax_length, ycm,
						t, form, labelpos, 0, fp );
			}
			}
		}
	} else {   /* log scale */
		if( dy > 0 ) {
		for( tref = ymin, t = ymin; t <= ymax; t += dy ) {
		psgetpos( x, t, &xcm, &ycm );
		psmove_cm( fp, xcm+ax_length, ycm );
		psdraw_cm( fp, xcm-ax_length, ycm );
		if( t > tref*10 ) {
			dy *= 10;
			tref = t;
		}
		}
		}
		fprintf( fp, "closepath stroke\n" );
		if(!( strcmp( which, "n" ) == 0  
			|| strcmp( which, "N") == 0 )) {
		if( d_label > 0 ) {
		for( t =label_min; t <= label_max; t *= d_label) {
			psgetpos(x, t, &xcm, &ycm ); 
			psnumberf_cm( xcm + ls * 1.3 * ax_length, ycm,
					t, form, labelpos, 0, fp );
			}
		}
		}
	}



	
	
}


void autoscale( float *x, int num,
	float *xmax, float *xmin, float *dx, float *Label_dx)
{
/*    scale x[0] --- x[num-1] ;
		max , min    2 digits max and min
		d            appropriate dx or dy  for axis
		Num_d                              for numbering
*/

	float max, min, d, Num_d;
	float tmax, tmin,  range;
	int i, t, dnum;
	float base,   frag;
	

	/* find real maximun and minimun		*/
	tmin = tmax = *x;
	for(i = 1; i < num; i++) {
		if( *(x+i) > tmax ) tmax = *(x+i);
		if( *(x+i) < tmin ) tmin = *(x+i);
	}


	range = tmax - tmin;
	if( range == 0 ) {
		max = min = tmax;
	} else {
		t = (int)log10( range ) ;
		base = pow( 10, t-1 );
	}


	if( tmax == 0 ) { max = 0;}
	else {
		if( tmax < 0 ) {	frag = -1;
					tmax = -tmax;
		} else {		frag = 1;
		}
		max  =  base * ceil( frag * tmax / base ) ;
	}



	if( tmin == 0 ) { min = 0;}
	else {
		if( tmin < 0 ) {	frag = -1;
					tmin = -tmin;
		} else {		frag = 1;
		}
		t = (int)log10( tmin );
		min  =  base * floor( frag * tmin / base ) ;
	}

	 if( max - min <= 0) { d = Num_d = 0; }
	 else {
		 t = (int)log10(max - min);
		 Num_d = d = pow( 10, t-1);
		 dnum = ( max - min )/ d ;
		 
		 if( dnum  > 70 )  {
		 	d *= 10; Num_d *=20;}
		 if( 70 >= dnum && dnum > 25 ) {
		 	d *= 5; Num_d *= 10;}
		 if( 25 >= dnum && dnum > 15 ) {
		 	Num_d *= 5;}
		 if( 15 >= dnum && dnum > 5 )  {
		 	Num_d *= 2;}

	}

	*xmax = max;  *xmin = min; *dx = d; *Label_dx = Num_d;
}


void autoscalelog( float *x, int num,
	float *xmax, float *xmin, float *dx, float *Label_dx)
{
/*    scale x[0] --- x[num-1] ;
		max , min    2 digits max and min
		d            appropriate dx or dy  for axis
		Num_d                              for numbering
		log version!
*/

	float max, min, d, Num_d;
	float tmax, tmin,  range;
	int i, t, dnum;
	float base,   frag;
	

	/* find real maximun and minimun		*/
	tmin = tmax = *x;
	for(i = 1; i < num; i++) {
		if( *(x+i) > tmax ) tmax = *(x+i);
		if( *(x+i) < tmin ) tmin = *(x+i);
	}

	if( tmin*tmax < 0 ) {
		fprintf( stderr, "different sign!!\n");
	}

	min = pow( 10, floor( log10(tmin) ) );
	max = pow( 10, ceil(  log10(tmax) ) );

	range = tmax - tmin;
	t = (int)log10( range ) ;
	if( max/min < 20 ) {
		base = pow( 10, t-1 );

		max  =  base * ceil( frag * tmax / base ) ;
		min =   10*base * floor( frag * tmax / (10*base) );
	}

	d = pow(10, floor( log10(tmin) ) );
	t = (int)log10(max/min);
	if( t < 1 ) {
		Num_d = 2;
	} else if( t < 4 ) {
		Num_d = 10;
	} else {
		Num_d = pow( 10, floor(t/4) );
	}


	*xmax = max;  *xmin = min; *dx = d; *Label_dx = Num_d;
}





void autoscale2( float *x, int num,
	float *xmax, float *xmin, float *dx, float *Label_dx, char fmt[])
{
/*    scale x[0] --- x[num-1] ;
		max , min    2 digits max and min
		d            appropriate dx or dy  for axis
		Num_d                              for numbering
*/

	float max, min, d, Num_d;
	float tmax, tmin,  range;
	int i, t, dnum;
	float base,   frag;
	float dmaxmin;
	

	/* find real maximun and minimun		*/
	tmin = tmax = *x;
	for(i = 1; i < num; i++) {
		if( *(x+i) > tmax ) tmax = *(x+i);
		if( *(x+i) < tmin ) tmin = *(x+i);
	}


	range = tmax - tmin;
	if( range == 0 ) {
		max = min = tmax;
	} else {
		t = (int)log10( range ) ;
		base = pow( 10, t-1 );
	}


	if( tmax == 0 ) { max = 0;}
	else {
		if( tmax < 0 ) {	frag = -1;
					tmax = -tmax;
		} else {		frag = 1;
		}
		max  =  base * ceil( frag * tmax / base ) ;
	}



	if( tmin == 0 ) { min = 0;}
	else {
		if( tmin < 0 ) {	frag = -1;
					tmin = -tmin;
		} else {		frag = 1;
		}
		t = (int)log10( tmin );
		min  =  base * floor( frag * tmin / base ) ;
	}

	 if( max - min <= 0) { d = Num_d = 0; }
	 else {
		 t = (int)log10(max - min);
		 Num_d = d = pow( 10, t-1);
		 dnum = ( max - min )/ d ;
		 
		 if( dnum  > 70 )  {
		 	d *= 10; Num_d *=20;}
		 if( 70 >= dnum && dnum > 25 ) {
		 	d *= 5; Num_d *= 10;}
		 if( 25 >= dnum && dnum > 15 ) {
		 	Num_d *= 5;}
		 if( 15 >= dnum && dnum > 5 )  {
		 	Num_d *= 2;}

	}

	*xmax = max;  *xmin = min; *dx = d; *Label_dx = Num_d;

	dmaxmin = max - min;
	if( dmaxmin <= 0.001 ) strcpy( fmt, "%.2e" );
	if( 0.001 < dmaxmin && dmaxmin <= 0.01 ) strcpy( fmt, "%.4f" );
	if( 0.01 < dmaxmin && dmaxmin <= 0.1 ) strcpy( fmt, "%.3f" );
	if( 0.1 < dmaxmin && dmaxmin <= 1 ) strcpy( fmt, "%.2f" );
	if( 1 < dmaxmin && dmaxmin <= 1000 ) strcpy( fmt, "%.1f" );
	if( 1000 < dmaxmin && dmaxmin  ) strcpy( fmt, "%.2e" );
}

	
void psputimage(int *image, int numx, int numy,
		float xmin, float xmax, float ymin, float ymax, FILE *fp)
{ /* put image  rgb * numx * numy 	image must be 3*numx*numy butes */


	int i, j, k;

	fprintf( fp, "/picstr %d string def\n", numx*3 );
	fprintf( fp, "%.2f %.2f translate\n", 
		psgraph.factor_xa * xmin + psgraph.factor_xb,
		psgraph.factor_ya * ymin + psgraph.factor_yb);
	
	fprintf( fp, "%.2f %.2f scale\n",
		psgraph.factor_xa*(xmax-xmin), psgraph.factor_ya*(ymax-ymin) );
	fprintf( fp, "%d %d 8\n", numx, numy);
	fprintf( fp, "[%d 0 0 %d 0 0]\n", numx, numy );	/*  ??? */
	fprintf( fp, "{currentfile picstr readhexstring pop}\n");
	fprintf( fp, "false 3 colorimage\n");
	for(i=0; i<numy; i++) {
		for(j=0; j<numx; j++) {
			for(k=0; k<3; k++) {
			if( *image > 15 ) {
				fprintf( fp, "%x", *image);
			} else {
				fprintf( fp, "0%x", *image);
			} 
			image++;
			}

		} fprintf( fp, "\n" );
	}

	fprintf( fp, "0 0 translate 1 1 scale\n");
}


int * data2rgb( float *data, int numx, int numy, float min, float max, int *rgb)
{
	/* data [ min: max ] ===>   blue -- green -- red */

	int i, j,k, *r, *g, *b;
	float d, f;

	r = rgb; 
	g = r + 1;
	b = r + 2;
	f = 2 / (max-min);

	for(i=0; i<numy; i++ ) {
		for(j=0; j<numx; j++ ) {
			d = f * (*(data+numx*i+j) - min);
			if( d < 1 ) {
				*r = 0;
				*g = rint(255 * d );
				*b = rint( 255 * ( 2 -2*d) );
			} else {
				*r = rint( 255 * 2 *(d-1) );
				*g = rint( 255 * (2-d) );
				*b = 0;
			}

		/*
			*g = rint( 255 * (-cos(d)+1) );
			if( d >M_PI ) {
				*r = rint( 255 * (-cos(d+M_PI)+1));
			} else {
				*r = 0;
			}
			if( d < M_PI ) {
				*b = rint( 255 * (cos(d)+1));
			} else {
				*b = 0;
			}
		*/
			r += 3; g +=3; b +=3;
		}
	}

	for(i=0; i<numy; i++ ) {
		for(j=0; j<numx; j++ ) {
			for(k=0; k<3; k++ ) {
				if( *(rgb+k+3*j+3*numx*i) < 0 )
					*(rgb+k+3*j+3*numx*i) = 0;
				if( *(rgb+k+3*j+3*numx*i) > 255 )
					*(rgb+k+3*j+3*numx*i) = 255;
				}
			}
		}


	return rgb;
}

int * data2rgb8( float *data, int numx, int numy, float min, float max, int *rgb)
{
	/* data [ min: max ] ===>   blue -- green -- red   8 level */

	int i, j,k, d, *rgbp;
	float f;
	int clo[8][3] = { 0, 0, 255,
			  135, 206, 235,
			  152, 251, 152,
			  0, 255, 0,
			  173, 255, 47,
			  255, 255, 0,
			  255, 165, 0,
			  255, 0, 0};
	rgbp = rgb;

	f = 8 / (max-min);

	for(i=0; i<numy; i++ ) {
		for(j=0; j<numx; j++ ) {
			d = (int)(f * (*(data+numx*i+j) - min));
			if( d > 8 )  d = 8;
			if( d < 0 )  d = 0;

			*rgb = clo[d][0]; rgb++;
			*rgb = clo[d][1]; rgb++;
			*rgb = clo[d][2]; rgb++;
		}
	}

	return rgbp;
}