ftp.nice.ch/pub/next/science/mathematics/hippoplotamus.2.0.s.tar.gz#/hippo2.0/hutil.c

This is hutil.c in view mode; [Download] [Up]

#include <string.h>
#include <stdlib.h>
#include <math.h>
#include <float.h>
#include "hippo.h"
#include "hutil.h"

GLOB_QUAL const char hippoutil_c_rcsid[] = "$Id: hutil.c,v 5.0 1993/08/17 21:56:48 rensing Exp $";

#define FLT_EQUAL(x,y) ((float)fabs((x)-(y))<=2.0*((y)*FLT_EPSILON+FLT_MIN))
#define MAX(x, y)  (((x) > (y)) ? (x) : (y))

/*
 * Internal functions
 */

static float calcTicks(float size, float *magnitude);

/*
 * Structures for the function registry.
 */

struct func_strt { 
     char name[FUNCNAMELEN+1]; 
     void *func;
};

static struct func_strt f_list1[] = {
{ "h_cut_lt", (void *)h_cut_lt },
{ "h_cut_gt", (void *)h_cut_gt },
{ "h_cut_le", (void *)h_cut_le },
{ "h_cut_ge", (void *)h_cut_ge },
{ "h_cut_inside", (void *)h_cut_inside },
{ "h_cut_outside", (void *)h_cut_outside },
{ "h_cut_in_incl", (void *)h_cut_in_incl },
{ "h_cut_out_incl", (void *)h_cut_out_incl },
{ "", (void *)NULL } };		/* must terminate with NULL */

static struct func_strt *f_list2 = NULL;

static int func_list_l = 0;


/*
 * Search for the function name, return its pointer.
 *   Return NULL if not found.
 */
void *h_fNameSrch(const char *fname)
{
     int i = 0;
     
     while( f_list1[i].func != NULL )
     {
	  if (strncmp(f_list1[i].name,fname,FUNCNAMELEN) == 0)
	       return f_list1[i].func;
	  i++;
     }

     if (f_list2 != NULL)
     {
	  i = 0;
	  while( f_list2[i].func != NULL )
	  {
	       if (strncmp(f_list2[i].name,fname,FUNCNAMELEN) == 0)
		    return f_list2[i].func;
	       i++;
	  }
     }
     
     return NULL;
}

/*
 * Search for the function pointer, return its name if found.
 */
char *h_fPtrSrch(void *func)
{
     int i = 0;
     
     while( f_list1[i].func != NULL )
     {
	  if (f_list1[i].func == func)
	       return f_list1[i].name;
	  i++;
     }
     
     if (f_list2 != NULL)
     {
	  i = 0;
	  while( f_list2[i].func != NULL )
	  {
	       if (f_list2[i].func == func)
		    return f_list2[i].name;
	       i++;
	  }
     }

     return NULL;
}


/*
 * Register the function on the registry.
 */
int h_func_reg(const char *fn, void *func)
{
     int i=0;
     char *sname;
     void *sfunc;
     
     if ((sname = h_fPtrSrch(func)) != NULL)
	  if (sname != fn)
	  {
	       h_error("h_funcReg: function name is in registry, but pointer does not match.");
	       return -1;
	  }
     	  else
	       return 0;

     if ((sfunc = h_fNameSrch(fn)) != NULL)
	  if (sfunc != func)
	  {
	       h_error("h_funcReg: function pointer is in registry, but name does not match.");
	       return -1;
	  }
     	  else
	       return 0;

     if (func_list_l == 0)
     {
	  f_list2 = (struct func_strt *)
	       malloc( (size_t)5*sizeof(struct func_strt) );
	  if (f_list2 == NULL)
	  {
	       h_error("h_funcReg: Error allocating memory for function registry.");
	       return -1;
	  }
	  func_list_l = 5;
	  i = 0;
     }
     else
     {
	  while( f_list2[i].func != NULL ) i++;
	  if (i==(func_list_l-1))
	  {
	       f_list2 = (struct func_strt *)
		    realloc( f_list2, 
			    (size_t)(func_list_l+5)*sizeof(struct func_strt) );
	       if (f_list2 == NULL)
	       {
		    h_error("h_funcReg: Error reallocating memory for function registry");
		    return -1;
	       }
	       func_list_l += 5;
	  }
     }

     strncpy(f_list2[i].name,fn,FUNCNAMELEN);
     f_list2[i].func = func;

     strncpy(f_list2[i+1].name,"",FUNCNAMELEN);
     f_list2[i+1].func = (void*)NULL;
     
     
     return 0;
}

void h_adjustAxis( float *low, float *high, int nbins, int log )
{
     int i;
     int axis = 0;
     float step, w, mag, range;
     float mylow, myhigh;
#define N_NICE 6
#ifndef __STDC__
     static
#endif
     float nice[N_NICE] = { 1.0,2.0,2.5,4.0,5.0,7.5 };
	  
     if (*low >= *high) 
     {
	  if (*low < 0.0) *high = 0.0;
	  else if (*low > 0.0) *low = 0.0;
	  else 
	  {
	       *low = -1.0;
	       *high = 1.0;
	       return;
	  }
     }
     
     if (nbins <= 0) 
     {
	  axis = 1;
	  nbins = 10;
     }
     
     /*
      * Round the "bin width" to a nice number.
      * If this is being done for an axis (ie nbins was 0), then
      * we don't have to go > *high.
      */
     w = (*high - *low)/((float)nbins);
     mag = floor(log10(w));
     i = 0;
     do
     {
 	  step = nice[i] * pow(10.0,mag);
     
	  mylow = floor(*low/step) * step;
	  myhigh = mylow + step*nbins;	  

	  i++;
	  if (i>=N_NICE) 
	  {
	       i = 0;
	       mag++;
	  }
     }
     while ( myhigh <= *high || FLT_EQUAL(*high,myhigh) );

     range = myhigh - mylow;
     if (axis && range < 1.1*(*high - *low)) range *= 1.1;
     
     /*
      * we now have decided on a range. Try to move low/high a little
      *  to end up on a nice number.
      *
      * first check if either end is near 0.0
      */
     if (!log && *low >= 0.0 && range>(1.05* *high) )
     {
	  *low = 0.0;
	  *high = range;
	  return;
     }
     if (*high<=0.0 && -range<(1.05* *low))
     {
	  *high = 0.0;
	  *low = -range;
	  return;
     }

     /*
      * try to round *low.
      */
     i = N_NICE-1;
     if (myhigh != 0.0)
	  mag = ceil(log10(fabs(myhigh)));
     else
	  mag = ceil(log10(fabs(mylow)));
     
     do
     {
	  step = nice[i] * pow(10.0,mag);
     
	  mylow = floor(*low/step) * step;
	  myhigh = mylow + range;

	  i--;
	  if (i<0) 
	  {
	       i = N_NICE-1;
	       mag--;
	  }
     } 
     while ((log && mylow <= 0.0) || myhigh <= *high );

     *low = mylow;
     *high = myhigh;
}
#undef N_NICE


/*
 * Hippo error message facility.
 */
void h_errmsg(const char *file, int line, const char *msg )
{
     fprintf(stderr, "Error in hippo file %s, line %d :\n", file, line);
     fprintf(stderr, "     %s\n", msg );
}


/*
 * Title expansion routine. Parse string for such things as 
 * %t, %x, %y, etc.
 */
char *h_expandLabel( char *dest, const char *src, int max, display disp )
{
     int i,j=0;
     int bind;
     int l;
     int k;
     float del,del2;
     char form[15];
     
     l = strlen(src);

     if (disp->tuple == NULL)
     {
	  strncpy( dest, src, max );
	  return dest;
     }
     
     for (i=0; i<l && j<max; i++)
     {
	  if (src[i] == '%')
	  {
	       i++;
	       switch (src[i])
	       {
	       case '%':
		    dest[j++] = '%';
		    break;
		    
	       case 't': case 'T':
		    if (disp->tuple->title != NULL)
		    {
			 strncpy(&(dest[j]),disp->tuple->title,max-j-1);
			 j = strlen(dest);
		    }
		    break;
		    
	       case 'x': case 'X':
		    if ((bind = disp->binding.x) >= 0)
		    {
			 strncpy(&(dest[j]),disp->tuple->label[bind],max-j-1);
			 j = strlen(dest);
		    }
		    break;

	       case 'y': case 'Y':
		    if ((bind = disp->binding.y) >= 0)
		    {
			 strncpy(&(dest[j]),disp->tuple->label[bind],max-j-1);
			 j = strlen(dest);
		    }
		    break;

	       case 'z': case 'Z':
		    if ((bind = disp->binding.z) >= 0)
		    {
			 strncpy(&(dest[j]),disp->tuple->label[bind],max-j-1);
			 j = strlen(dest);
		    }
		    break;
		    
	       case 'w': case 'W':
		    if ((bind = disp->binding.weight) >= 0)
		    {
			 strncpy(&(dest[j]),disp->tuple->label[bind],max-j-1);
			 j = strlen(dest);
		    }
		    break;
		    
	       case 'd': case 'D':
		    i++;
		    switch (src[i])
		    {
		    case 'x': case 'X':
			 del = (disp->xAxis.high-disp->xAxis.low)
			      /(float)disp->bins.xAxis.nBins;
			 for (k=1; k<6; k++)
			 {
			      del2 = del*pow(10.0,(double)k);
			      if (fabs((int)del2 - del2) < 1e-2) break;
			 }
			 sprintf(form,"%%1.%df",k);
			 sprintf(&dest[j],form,del);
			 j = strlen(dest);
			 break;
			 
		    case 'y': case 'Y':
			 del = (disp->yAxis.high-disp->yAxis.low)
			      /(float)disp->bins.yAxis.nBins;
			 for (k=1; k<6; k++)
			 {
			      del2 = del*pow(10.0,(double)k);
			      if (fabs((int)del2 - del2) < 1e-2) break;
			 }
			 sprintf(form,"%%1.%df",k);
			 sprintf(&dest[j],form,del);
			 j = strlen(dest);
			 break;
		    }
		    break;
		    
	       case 'e': case 'E':
		    i++;
		    switch (src[i])
		    {
		    case 'x': case 'X':
			 if ((bind = disp->binding.xerror) >= 0)
			 {
			      strncpy(&(dest[j]),disp->tuple->label[bind],
				      max-j-1);
			      j = strlen(dest);
			 }
			 break;
			 
		    case 'y': case 'Y':
			 if ((bind = disp->binding.yerror) >= 0)
			 {
			      strncpy(&(dest[j]),disp->tuple->label[bind],
				      max-j-1);
			      j = strlen(dest);
			 }
			 break;
		    }
		    break;
	       }
	  }
	  else
	       dest[j++] = src[i];
     }

     if (j<max)
	  dest[j] = '\0';
     else 
	  dest[max-1] = '\0';

     return dest;
}


/*
 * Standard Cut functions.
 */

int h_cut_lt( float **ntdata, int row, double param[] )
{
     if (ntdata[(int) param[0] ][row] < param[1])
	  return 1;
     else
	  return 0;
}

int h_cut_gt( float **ntdata, int row, double param[] )
{
     if (ntdata[(int) param[0]][row] > param[1] )
	  return 1;
     else
	  return 0;
}

int h_cut_le( float **ntdata, int row, double param[] )
{
     if (ntdata[(int) param[0]][row] <= param[1] )
	  return 1;
     else
	  return 0;
}

int h_cut_ge( float **ntdata, int row, double param[] )
{
     if (ntdata[(int) param[0]][row] >= param[1] )
	  return 1;
     else
	  return 0;
}

int h_cut_inside( float **ntdata, int row, double param[] )
{
     if (ntdata[(int) param[0]][row] > param[1] &&
	 ntdata[(int) param[0]][row] < param[2])
	  return 1;
     else
	  return 0;
}

int h_cut_outside( float **ntdata, int row, double param[] )
{
     if (ntdata[(int) param[0]][row] < param[1] ||
	 ntdata[(int) param[0]][row] > param[2])
	  return 1;
     else
	  return 0;
}

int h_cut_in_incl( float **ntdata, int row, double param[] )
{
     if (ntdata[(int) param[0]][row] >= param[1] &&
	 ntdata[(int) param[0]][row] <= param[2])
	  return 1;
     else
	  return 0;
}

int h_cut_out_incl( float **ntdata, int row, double param[] )
{
     if (ntdata[(int) param[0]][row] <= param[1] ||
	 ntdata[(int) param[0]][row] >= param[2])
	  return 1;
     else
	  return 0;
}


int doCuts( float **ntdata, int row, func_id cutlist )
{
     typedef int (*cutfunc)(float **, int, double *);
	  
     func_id c = cutlist;
     int result = 0;
     
     while (c != NULL)
     {
	  result += (*((cutfunc)c->funcPtr))(ntdata, row,
					     (double *)c->paramBlk );
	  c = c->next;
     }

     return result;
}


int genTicks(axis_t *thisaxis, int maxticks, int *nticks, 
	float *ticks, char labels[][10], float *pmag) 
{
     int i;
     
     float magnitude,y=0.0,yr,startTick,tickSize;
     char string[70];
     char pstr[10];
     int nLogTicks;
     float logTicks[5];
     float maglow, maghigh, magrng, magStep;
     
     *nticks = 0;
     
     if (thisaxis->low >= thisaxis->high)
     {
	  h_error("h_plot: axis low > high");
	  sprintf(string,"low = %f, high = %f\n",thisaxis->low,
		  thisaxis->high);
	  h_error(string);
	  return -1;
     }

     if (!thisaxis->flags.log)
     {
	  tickSize = calcTicks((thisaxis->high -
				thisaxis->low), &magnitude);
	  startTick = ceil( thisaxis->low / tickSize) * tickSize;

	  if (fabs(magnitude) <= 3.0)
	       *pmag = 0.0;
	  else 
	  {
	       if (startTick != 0.0)
		    *pmag = floor(log10(fabs(startTick)));
	       else
		    *pmag = magnitude;
	  }
	  	       	  
	  sprintf(pstr,"%%1.%df",(int)MAX((*pmag-magnitude),0.0));

	  y = startTick;
	  while (y <= thisaxis->high || FLT_EQUAL(thisaxis->high,y))
	  {
	       if (*nticks >= maxticks)
	       {
		    h_error("Too many ticks along y axis.");
		    return -1;
	       }
	       
	       yr = floor(y/pow(10.0,magnitude) + 0.5);
	       ticks[*nticks] = yr * pow(10.0,magnitude);
	       sprintf(labels[*nticks],pstr,yr*pow(10.0,magnitude-*pmag));
	       (*nticks)++;
	       y += tickSize;
	  }
	  if (fabs(magnitude) <= 3.0) magnitude = 0.0;
     }
     else
     {
	  maghigh = ceil(log10(thisaxis->high)*(1.0-2.0*FLT_EPSILON));
	  maglow = floor(log10(thisaxis->low)*(1.0+2.0*FLT_EPSILON));
	  magrng = log10(thisaxis->high/thisaxis->low);
	  
	  if (magrng <4.0)
	  {
	       nLogTicks = 3;
	       logTicks[0] = 1.0;
	       logTicks[1] = 2.0;
	       logTicks[2] = 5.0;
	       magStep = 1.0;
	  }
	  else if (magrng < 7.0)
	  {
	       nLogTicks = 2;
	       logTicks[0] = 1.0;
	       logTicks[1] = 2.0;
	       magStep = 1.0;
	  }
	  else
	  {
	       nLogTicks = 1;
	       logTicks[0] = 1.0;
	       magStep = 2.0;
	  }
	       
	  if (nLogTicks > 1 && (fabs(maglow)>3 || fabs(maghigh)>3))
	       *pmag = maglow;
	  else
	       *pmag = 0.0;
	  magnitude = maglow;
	  i = 0;
	  while ((y=logTicks[i]*pow(10.0,magnitude)) <= thisaxis->high 
		 || FLT_EQUAL(thisaxis->high,y))
	  {
	       if (*nticks >= maxticks)
	       {
		    h_error("Too many ticks along y axis.");
		    return -1;
	       }
	       if (y >= thisaxis->low)
	       {
		    ticks[*nticks] = log10(y);
		    /*
		     * be careful: there is a bug in the NeXT (s)printf 
		     *   routine when you do, eg. printf("%1.0g",0.01);
		     */
		    if ((magnitude-*pmag) > 4 || (magnitude-*pmag) < -3)
			 strcpy(pstr,"%1.0e");
		    else
			 sprintf(pstr,"%%1.%df",
				 (int)((magnitude-*pmag)>0?
				       0:-(magnitude-*pmag)));
		    sprintf(labels[*nticks],pstr,y*pow(10.0,-*pmag));
		    (*nticks)++;
	       }
	       
	       i++;
	       if (i>=nLogTicks)
	       {
		    i = 0;
		    magnitude += magStep;
	       }
	  }
    }
    return 0;
}     

#define MIN_TICKS 4

static float calcTicks(float size, float *magnitude)
{
     static float	goodTicks[] = {10.0, 5.0, 4.0, 2.0, 1.0};
     float		tickSize;
     int		tickIndex;

     if (size <= 0.0) 
     {
	  h_error("calcTicks: size is <= 0");
	  size = fabs(size);
	  if (size == 0.0) size = 1.0;
     }
     
     *magnitude = floor(log10(size));
     if (size/pow(10.0,*magnitude) < MIN_TICKS) (*magnitude)--;

     /*
      *  now fit the max number of ticks into this range
      */
     for (tickIndex = 0;
	  size / (tickSize=goodTicks[tickIndex]*pow(10.0,*magnitude) ) 
	  < MIN_TICKS;
	  tickIndex++);

     if (tickIndex == 0) 
	  (*magnitude)++;
     
     return tickSize;
}

These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.