Return to capaUnit.c CVS log

Up to [LON-CAPA] / capa / capa51 / pProj

- BUG#3952, the instructor's unit was being parsed by the stupid, but fast unit parser, swapping it out so that the instructor's parser is the same as the one used to parse the students's answer

/* functions to handle the unit parser/comparison engine
   Copyright (C) 1992-2000 Michigan State University

   The CAPA system is free software; you can redistribute it and/or
   modify it under the terms of the GNU General Public License as
   published by the Free Software Foundation; either version 2 of the
   License, or (at your option) any later version.

   The CAPA system is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
   General Public License for more details.

   You should have received a copy of the GNU General Public
   License along with the CAPA system; see the file COPYING.  If not,
   write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
   Boston, MA 02111-1307, USA.

   As a special exception, you have permission to link this program
   with the TtH/TtM library and distribute executables, as long as you
   follow the requirements of the GNU GPL in regard to all of the
   software in the executable aside from TtH/TtM.
*/

/* =||>|===================== capaUnit.c   =====================|<||= */
/*   created by Isaac Tsai   1997                                    */
/*   by Isaac Tsai 1997, 1998, 1999                      */
/* =||>|========================================================|<||= */
#include <stdio.h>        /* fopen()  */
#include <stdlib.h>
#include <ctype.h>        /* isalnum()   */
#include <string.h>
#include <math.h>
#include <float.h>

#include "capaParser.h"

int      PrefixTbl[QUARTER_K];
int      BaseUnitcnt;
double   CScale[BASEUNIT_LIMIT];
double   CExp[BASEUNIT_LIMIT];
char     CSymb[BASEUNIT_LIMIT][SYMBOL_MAXLEN];
Unit_t  *UnitTree_p;
double   MinSquared;
Unit_t  *MinSquaredUnit_p;
Unit_t  *InqueryUnit_p;
double  *TmpAexp, *TmpBexp;
Unit_t  *EquivUnit[BASEUNIT_LIMIT];
double   MinValue[BASEUNIT_LIMIT];
int      EquivUnitCnt;
char     Sbuf[ONE_K_SIZE];
int      Sidx;
Unit_t  *Pstack[ONE_K_SIZE];
int      Ptopidx;
int      gUnitError;

FILE    *ufp;

/* ==================================================================== */
void c_ignorewhite(FILE *f) /* ignore white spaces from a file stream */
{
  register int c;
  register int ok;
 
  ok = 0;
  do {
    do {  c = getc(f);
    } while ( isspace(c) );
    ungetc(c,f);
    if (c == '#') {
      while (getc(f) != '\n');
    } else ok = 1;
  } while( ! ok);
}

int c_getint(FILE *f)  /* returns an integer from the file stream */
{
    int c;
    int value;
 
    c_ignorewhite(f);
    c = fgetc(f);
    if (!isdigit(c)) {
        fprintf(stderr,"Error: Expected digit, got %c\n", c);
        exit(-1);
    }
    ungetc(c,f);
    fscanf(f,"%d", &value);
    return(value);
}
int c_getsec_range(FILE *f,int *low,int *high)
{
    int c;
    int tmp, result;
 
    c_ignorewhite(f);
    c = fgetc(f);
    if( c == '[' ) { /* specify a range of sections */
      do {  c = getc(f); } while ( isspace(c) );
      if (!isdigit(c)) {
        fprintf(stderr,"Error in section range format, expecting a number.\n");
        result = -1;
        return (result);
      }
      ungetc(c,f);
      fscanf(f,"%d", low);
      do {  c = getc(f); } while ( isspace(c) );
      if( c == ',' ) {
        do {  c = getc(f); } while ( isspace(c) );
        if (!isdigit(c)) {
          fprintf(stderr,"Error in section range format, expecting a number.\n");
          result = -1;
          return (result);
        }
        ungetc(c,f);
        fscanf(f,"%d", high);
        do {  c = getc(f); } while ( isspace(c) );
        if( c == ']' ) {
          if( *high < *low ) {
            tmp= *high; *high = *low; *low =tmp;
          }
          if(*low <=0) {
            *low = 1;
          }
          if(*high <=0) {
            *high =1;
          }
          /* printf("Section range=>[%d,%d]\n",*low,*high); */
          result = 2;
        }
      } else { /* no , specified */
        result = -1;
        return (result);
      }
    } else { /* specify a section only */
      if (!isdigit(c)) {
        fprintf(stderr,"Error: Expected digit, got %c\n", c);
        result = -1;
        return (result);
      }
      ungetc(c,f);
      fscanf(f,"%d", low);
      result = 1;
    }
    return (result);
}

double c_getdouble(FILE *f)
{
    int c;
    double value;
 
    c_ignorewhite(f);
    c = fgetc(f);
    if (!isdigit(c)) {
        fprintf(stderr,"Error: Expected digit, got %c\n", c);
        exit(-1);
    }
    ungetc(c,f);
    fscanf(f,"%lf", &value);
    return(value);
}

/*   read until encountered an unrecognizable char */
/*      space, #, anything other than alphanum, {}-^_ */
char *c_getword(FILE *f) 
{
  register int c;
  register int idx;
  char     tmp_string[ONE_K];
  char     *new_string;

  idx = 0;
  c_ignorewhite(f);
    do {  c = getc(f);
      tmp_string[idx] = c;
      idx++;
    } while (isalnum(c) || c == '{' || c == '}' || c == '-' || 
             c == '^'   || c == '_' );
    ungetc(c,f); idx--;
    tmp_string[idx] = 0;
    new_string = (char *)malloc( (idx+1)*sizeof(char) );
    strncpy(new_string,tmp_string, (idx+1) );
  
  return (new_string);
}
/*   read until encountered a newline, # */
char *c_getstring(FILE *f) 
{
  register int c;
  register int idx;
  char     tmp_string[1024];
  char     *new_string;

  idx = 0;
  c_ignorewhite(f);
    do {  c = getc(f);
      tmp_string[idx] = c;
      idx++;
    } while (isalnum(c) || c == '{' || c == '}' || c == '-' || 
             c == '^'   || c == ' ' || c == ',' || c == ';' ||
             c == '.'   || c == '(' || c == ')' || c == '=' ||
             c == '+'   || c == '*' || c == '/' );
    ungetc(c,f); idx--;
    tmp_string[idx] = 0;
    c = tmp_string[idx-1];
    while( c == ' ') {   /* get rid of trailing white space */
       idx--;
       c = tmp_string[idx-1];
    }
    tmp_string[idx] = 0;
    new_string = (char *)malloc( (idx+1)*sizeof(char) );
    strncpy(new_string,tmp_string, (idx+1) );
  
  return (new_string);
}
char *c_getcomment(FILE *f) 
{
  register int c;
  register int idx;
  char     tmp_string[ONE_K];
  char     *new_string;

  idx = 0;
  while (getc(f) != '#');
  while ((c = getc(f)) == ' ');  ungetc(c,f);
    do {  c = getc(f);
      tmp_string[idx] = c;
      idx++;
    } while ( isprint(c) );
/*
    } while (isalnum(c) || c == '{' || c == '}' || c == '-' || 
             c == '^'   || c == ' ' || c == ',' || c == ';' ||
             c == '.'   || c == '(' || c == ')' || c == '=' );
*/
    ungetc(c,f); idx--;
    tmp_string[idx] = 0;
    c = tmp_string[idx-1];
    while( c == ' ') {   /* get rid of trailing white space */
       idx--;
       c = tmp_string[idx-1];
    }
    tmp_string[idx] = 0;
    new_string = (char *)malloc( (idx+1)*sizeof(char) );
    strncpy(new_string,tmp_string, (idx+1) );
  
  return (new_string);
}
void  c_moveto_unit(FILE *f)
{
  register int c;
  register int ok;
 
  ok = 0;
  do {
    do {  c = getc(f);
    } while (c != '<' );
    c = getc(f);
    if (c == '<') {
      ungetc(c,f); ungetc(c,f); ok=1;
    }
  } while( ! ok);
}

int  c_gettype(FILE *f)
{
  register int c;
  register int idx;
  char     tmp_string[ONE_K];
  char     new_string[ONE_K];
  
  idx = 0;
PRESTART:
  c_ignorewhite(f);
  while ((c=getc(f)) != '<') { if ( (char)c==(char)EOF ) return U_UNKNOWN; }
  c = getc(f);
  if( c == '<' ) {
    c_ignorewhite(f);
PREEND:
    do {  c = getc(f);
      tmp_string[idx] = toupper(c);
      idx++;
    } while ( c != '>' );
    c = getc(f); 
    if( c == '>' ) {
      idx--;
      tmp_string[idx] = 0;
      c = tmp_string[idx-1];
      while( c == ' ') {   /* get rid of trailing white space */
        idx--;
        c = tmp_string[idx-1];
      }
      tmp_string[idx] = 0;
      strncpy(new_string,tmp_string, (idx+1) );
    } else {
      ungetc(c,f);
      goto PREEND;
    }
  } else {
    goto PRESTART;
  }
  if( !strcmp(new_string,"BASE UNIT") ) {
    return (U_BASE);
  }  
  if( strcmp(new_string, "DERIVED UNIT") == 0 ) {
    return (U_DERIVED);
  }
  if( strcmp(new_string, "PREFIX") == 0 ) {
    return (U_PREFIX);
  }
  if( strcmp(new_string, "CONSTANTS") == 0 ) {
    return (U_CONSTANT);
  }
  if( strcasecmp(new_string, "DEFAULTS") == 0 ) {
    return (U_DEFAULT);
  }
  return (U_UNKNOWN);
  
}

/* =================================================================== */
/* =================================================================== */
/* returns: 0 success */
/*          1 the first units string u1_str could not be reduce to a valid unit */
/*          2 the second units string could not be reduced to a valid unit */
int
u_convert_unit(char *u1_str,char *u2_str,double *ratio)
{
  Unit_t   *ap, *bp;
  int       result=0;
  
  while( isspace(*u1_str) )  u1_str++;
  while( isspace(*u2_str) )  u2_str++;
  bp = parse_unit_expr(u2_str);
  Ptopidx=0;
  postwalk_utree(bp);
  if( Ptopidx == 1 ) {
    simplify_unit(Pstack[Ptopidx]);
    bp = Pstack[Ptopidx];
    /* print_unit_t(bp); */
    ap = parse_unit_expr(u1_str);
    Ptopidx=0;
    postwalk_utree(ap);
    if( Ptopidx == 1 ) {
      simplify_unit(Pstack[Ptopidx]);
      /* print_unit_t(Pstack[Ptopidx]); */
      if( (Pstack[Ptopidx]->u_count != 0) ||
          (Pstack[Ptopidx]->u_count == bp->u_count) ) { /* has unit */
        *ratio = units_ratio(Pstack[Ptopidx], bp);
      } else {
        result = 1;
      }
    }
    free_utree(ap);
  } else {
    result = 2;
  }
  free_utree(bp);
  return (result);
}

/* =================================================================== */



Unit_t *
u_find_symb (char *name, Unit_t *t, int *result) 
{
  
  if (t == NULL)  return t;

  for (;;) {
    if ( comp_unit_symb(name,t->u_symbol) < 0 ) {
      if (t->u_left == NULL)  {
        /* printf("L not found\n"); */
        *result = 0;
        break;
      }
      t = t->u_left;
    } else if ( comp_unit_symb(name,t->u_symbol) > 0 ) {
      if (t->u_right == NULL) {
        /* printf("R not found\n"); */
        *result = 0;
        break;
      }
      t = t->u_right;
    } else {
     *result = 1;
      break;
    }
  }
  return t;
}
/* ------------------------------------------------------------- */
/*   use the input unit_t's element list to locate the min squared */
/*   error fit of the unit tree        */
/*   report either exact fit or approx */

void
u_find_name(Unit_t *t)
{
  int      ii;
  Unit_E  *eu_p;
  
  MinSquared = FLT_MAX;
  EquivUnitCnt=0;
  InqueryUnit_p = t;
  /* printf("INQ[[%s,%s,%d]]\n", U_SYMB(t), U_NAME(t),U_COUNT(t)); */
  TmpAexp = (double *)capa_malloc(BaseUnitcnt,sizeof(double));
  TmpBexp = (double *)capa_malloc(BaseUnitcnt,sizeof(double));
  for(ii=0;ii<BaseUnitcnt;ii++) {
     TmpAexp[ii] = 0.0;
  }
  if( t->u_count > 0 ) {
    for(eu_p = t->u_list; eu_p; eu_p = eu_p->ue_nextp) {
      TmpAexp[eu_p->ue_index] = eu_p->ue_exp;
      /* printf("(%d)^(%g) ",eu_p->ue_index,TmpAexp[eu_p->ue_exp]); */
    }
    /* printf("\n"); */
  }
  inorder_diff(UnitTree_p);
  /*capa_mfree((char *)TmpAexp); capa_mfree((char *)TmpBexp);*/
  
}

void
print_matches(Unit_t *t)
{
  double   scale, factor;
  Unit_t  *tmp_p;
  int      ii;
  
  scale = t->u_scale;
  if( MinSquared == 0.0 ) {  /* exact match */
    if( EquivUnitCnt > 0 ) {
      printf(" Entered unit is equivalent to:\n");
      for(ii=0;ii<EquivUnitCnt;ii++) {
        tmp_p = EquivUnit[ii];
        if( MinSquared ==  MinValue[ii] ) {
          if( tmp_p->u_type == U_BASE ) {    /* if there is a base unit */
            MinSquaredUnit_p = tmp_p;
          }
          factor = scale / tmp_p->u_scale;
          printf(" <<%g %s>>", factor,U_SYMB(tmp_p));
        }
      }
      printf("\n");
      
    }
  } else {  /* no exact match */
    if( EquivUnitCnt > 0 ) {
      printf(" Entered unit is approximated by:\n");
      for(ii=0;ii<EquivUnitCnt;ii++) {
        tmp_p = EquivUnit[ii];
        if( MinSquared ==  MinValue[ii] ) {
          printf(" <<%s>> ", U_SYMB(tmp_p) );
        }
      }
      printf("\n");
    }
  }
}

/* ------------------------------------ */
double
u_squared_diff(Unit_t  *a, Unit_t *b)
{
  double   result;
  double   squared_diff = 0.0;
  int      ii;
  Unit_E  *eu_p;
  
  
  for(ii=0;ii<BaseUnitcnt;ii++) {
    TmpAexp[ii] = 0.0;
    TmpBexp[ii] = 0.0;
  }
  if( a->u_count > 0 ) {
    for(eu_p= a->u_list; eu_p; eu_p = eu_p->ue_nextp) {
      TmpAexp[eu_p->ue_index] = eu_p->ue_exp;
    }
  }
  if( b->u_count > 0 ) {
    for(eu_p= b->u_list; eu_p; eu_p = eu_p->ue_nextp) {
      TmpBexp[eu_p->ue_index] = eu_p->ue_exp;
      /* printf("Exp[%d]=%g ",ii,TmpBexp[ii]); */
    }
    /* printf("\n"); */
  }
  for(ii=0;ii<BaseUnitcnt;ii++) {
    result = TmpAexp[ii] - TmpBexp[ii];
    squared_diff = squared_diff + result*result;
  }
  
  return (squared_diff);
}

double
u_sq_diff(Unit_t *b)
{
  double   result;
  double   squared_diff = 0.0;
  int      ii;
  Unit_E  *eu_p;
  
  
  for(ii=0;ii<BaseUnitcnt;ii++) {
    TmpBexp[ii] = 0.0;
  }
  if( b->u_count > 0 ) {
    for(eu_p= b->u_list; eu_p; eu_p = eu_p->ue_nextp) {
      TmpBexp[eu_p->ue_index] = eu_p->ue_exp;
      /* printf("Exp[%d]=%g ",ii,TmpBexp[ii]); */
    }
    /* printf("\n"); */
  } else if( b->u_type == U_BASE  ) {
    TmpBexp[b->u_index] = 1.0;
  }
  for(ii=0;ii<BaseUnitcnt;ii++) {
    result = TmpAexp[ii] - TmpBexp[ii];
    squared_diff = squared_diff + result*result;
  }
  
  return (squared_diff);

}
/* ------------------------------------ */

int
inorder_diff(node_p) Unit_t  *node_p;
{
  int      result;
  double   sq_diff=0.0;
  
  if( node_p == NULL )  return (1);
  
  result = inorder_diff(U_LEFT(node_p));
  if( result ) {
     sq_diff = u_sq_diff(node_p);
     /*
     printf("DIFF [%s,%s,%d] - [%s,%s,%d] = %g\n", 
      U_SYMB(InqueryUnit_p), U_NAME(InqueryUnit_p),U_COUNT(InqueryUnit_p),
      U_SYMB(node_p), U_NAME(node_p),U_COUNT(node_p),sq_diff);
     */
     if( MinSquared > sq_diff) {
       MinSquaredUnit_p = node_p;
       MinSquared = sq_diff;
     } else if ( MinSquared == sq_diff) {
       EquivUnit[EquivUnitCnt] = node_p;
       MinValue[EquivUnitCnt] = sq_diff;
       EquivUnitCnt++;
     }
  }
  result = inorder_diff(U_RIGHT(node_p));
  
  return (result);
}


int
alphaorder_utree(node_p) Unit_t  *node_p;
{
  int  result;
  
  if( node_p == NULL )  return (1);
  
  result = alphaorder_utree(U_LEFT(node_p));
  if( result ) printf(" (%s,%s)\n", U_SYMB(node_p), U_NAME(node_p) );
  result = alphaorder_utree(U_RIGHT(node_p));
  
  return (result);
}

int
w_alphaorder_utree(node_p) Unit_t  *node_p;
{
  int  result;
  
  if( node_p == NULL )  return (1);
  
  result = alphaorder_utree(U_LEFT(node_p));
  if( result ) { 
     printf(" (%s,%s)\n", U_SYMB(node_p), U_NAME(node_p) );
  }
  result = alphaorder_utree(U_RIGHT(node_p));
  
  return (result);
}

/* --------------------------------------------------------------------- */
void
print_unit_tree(int mode)
{
  if( mode == 1 ) {
    alphaorder_utree(UnitTree_p);
  } else {
    w_alphaorder_utree(UnitTree_p);
  }
}


int
preorder_utree(node_p) Unit_t  *node_p;
{
  int  result;
  
  if( node_p == NULL )  return (1);
  printf("Preorder=[[%s,%s,%d]]\n", U_SYMB(node_p), U_NAME(node_p),U_COUNT(node_p));
  result = preorder_utree(U_LEFT(node_p));
  if( result ) result = preorder_utree(U_RIGHT(node_p));
  return (result);
}
int
inorder_utree(node_p) Unit_t  *node_p;
{
  int  result;
  
  if( node_p == NULL )  return (1);
  
  result = inorder_utree(U_LEFT(node_p));
  if( result ) printf("INorder=[[%s,%s,%d]]\n", 
    U_SYMB(node_p), U_NAME(node_p),U_COUNT(node_p));
  result = inorder_utree(U_RIGHT(node_p));
  
  return (result);
}
int
postorder_utree(node_p) Unit_t  *node_p;
{
  int  result;
  
  if( node_p == NULL )  return (1);
  
  result = postorder_utree(U_LEFT(node_p));
  if( result ) result = postorder_utree(U_RIGHT(node_p));
  if( result ) {
    switch(U_TYPE(node_p)) {
      case U_DERIVED:   print_unit_t(node_p);
            break;
      case U_CONSTANT:  printf("(%g)",U_SCALE(node_p));
            break;
      case U_OP_POWER:  printf("^");
            break;
      case U_OP_TIMES:  printf("*");
            break;
      case U_OP_PLUS:   printf("+");
            break;
      case U_OP_MINUS:  printf("-");
            break;
      case U_OP_DIVIDE: printf("/");
            break;
      default:          printf("()");
            break;  
    }
  }
  return (result);
}

/* returns 1 on okay, 2 on error*/
int
postwalk_utree(Unit_t  *n_p)
{
  int  result;
  
  if( n_p == NULL )  return (1);
  
  result = postwalk_utree(U_LEFT(n_p));
  if (result !=2) {
    if( result ) result = postwalk_utree(U_RIGHT(n_p));
    if (result !=2) {
      if( result ) {
	switch(U_TYPE(n_p)) {
	case U_DERIVED:   Ptopidx++; Pstack[Ptopidx] = n_p;  /* push into stack */
	  break;
	case U_CONSTANT:  Ptopidx++; Pstack[Ptopidx] = n_p;  /* push into stack */
	  break;
	case U_UNKNOWN:   result=2; 
	  /*push into stack anyway, try to parse rest of tree */
	  break;
	case U_OP_POWER:  printf("^"); result=2;
	  break;
	case U_OP_TIMES:  process_op(U_OP_TIMES);        /* process operator */
	  break;
	case U_OP_PLUS:   printf("+"); result=2;
	  break;
	case U_OP_MINUS:  printf("-"); result=2;
	  break;
	case U_OP_DIVIDE: process_op(U_OP_DIVIDE);       /* process operator */
	  break;
	default:          printf("()"); result=2;
	  break;  
	}
      }
    }
  }
  return (result);
}

void
process_op(int op)
{
  Unit_t  *ap, *bp;
  double   exp_scale;
  int      no_error=1;
  
  bp = Pstack[Ptopidx--]; 
  ap = Pstack[Ptopidx--]; 
  
  switch(op) {
    case U_OP_TIMES:  exp_scale = 1.0;  break;
    case U_OP_DIVIDE: exp_scale = -1.0; break;
    case U_OP_PLUS:   
    case U_OP_MINUS:  no_error = u_pm_op(ap,bp,op);
                      if(no_error) {
                        Ptopidx++;
                        Pstack[Ptopidx] = ap;
                      }
                      break;
    default:          no_error=0; 
                      printf("No such op on the parse tree!\n");
          break;
  }
  if(no_error) {
    u_copy_unit(ap, bp, exp_scale);
    Ptopidx++;
    Pstack[Ptopidx] = ap;
  }
}

Unit_t*
process_utree(Unit_t *t)
{
  Ptopidx=0;
  postwalk_utree(t);
  if( Ptopidx == 1 ) {
    //fprintf(stderr,"Correctly parsed!\n");
    //fprintf(stderr,"Unit:%s\n",Sbuf);
    simplify_unit(Pstack[Ptopidx]);
    //Pstack[Ptopidx]->u_symbol[0]='\0';
    //fprintf(stderr,Pstack[Ptopidx]->u_symbol,"");
    print_unit_t(Pstack[Ptopidx]);
    //u_find_name(Pstack[Ptopidx]);
    //print_matches(Pstack[Ptopidx]);
    return(Pstack[Ptopidx]);
    //free_utree(t);
  }
  return(t);
}

/* ============================================================== */
/*  called from capaCommon.c */
/*                      */
/*  UNIT_FAIL           */
/*  NO_UNIT             */
/*  result: UNIT_OK correct   */
/*                            */
/* -------------------------------------------------------------- */
int  check_correct_unit(char *u_symb,Unit_t *t,double *scale)
{
  Unit_t   *ap;
  int       result=UNIT_OK;

#ifdef UNIT_DBUG
   if ((ufp=fopen("unit.DBUG","a"))==NULL) { fprintf(stderr,"Error: can't open login debug\n"); return UNIT_FAIL; }
#endif 

  while( isspace(*u_symb) )  u_symb++; 
  /* <= change this to search from the end of string */
  /* or to get rid of all the white spaces */


  ap = parse_unit_expr(u_symb);
  Ptopidx=0;

  if (postwalk_utree(ap)==1) {
#ifdef UNIT_DBUG
    fprintf(ufp,"Ptopidx %d\n",Ptopidx);
#endif
    if( Ptopidx == 1 ) {
      simplify_unit(Pstack[Ptopidx]);
      
      if( (Pstack[Ptopidx]->u_count != 0) ||
	  (Pstack[Ptopidx]->u_count == t->u_count) ) { /* has unit */
	*scale = units_ratio(Pstack[Ptopidx], t);
	if( *scale == 0.0 ) {
	  result = UNIT_IRRECONCIBLE;
	}
	free_utree(ap);
      } else {
	result = UNIT_INVALID_STUDENT3;
      }
    } else { /* invalid unit representation */
      result = UNIT_INVALID_STUDENT2;
    }
  } else {
    result = UNIT_INVALID_STUDENT1;
  }
#ifdef UNIT_DBUG
  fclose(ufp);
#endif 
  return (result);
}

/* ============================================================= */
int
free_units()
{
  free_utree(UnitTree_p);
  UnitTree_p=NULL;
  return 0;
}

int
free_utree(Unit_t  *t)
{
  int  result=1;
  
  if( t == NULL )  return (1);
  u_postfree(t);
  t=NULL;
  
  return (result);
}


int
u_postfree(Unit_t  *t)
{
  int  result;
  
  if( t == NULL )  return (1);
  
  result = u_postfree(U_LEFT(t));
  if( result ) result = u_postfree(U_RIGHT(t));
  if( result ) {
    if( t->u_comment ) {
      capa_mfree((char *)t->u_comment);
    }
    freelist_unit_e(t->u_list);
    capa_mfree((char *)t);
  }
  return (result);
}


void
print_unit_t(Unit_t *t) 
{
  Unit_E  *ue_p;

  /* printf("  Unit::[%s,%d]= %g * ", t->u_symbol,t->u_count,t->u_scale); */
  printf("  Unit::[%s] = %g * ", t->u_symbol, t->u_scale);
  for(ue_p=t->u_list; ue_p ; ue_p = ue_p->ue_nextp) {
    /*
    printf("<%s,%d,%g,%g> ",ue_p->ue_symbol,ue_p->ue_index,ue_p->ue_scale,ue_p->ue_exp);
    */
    printf("(%g*%s^%g) ",ue_p->ue_scale,ue_p->ue_symbol,ue_p->ue_exp);
  }
  printf("\n");

}
/*  ----------------------------------------------------------- */
/*  copy the Unit_E linked list from b_p->u_list to a_p->u_list */
/*   create some Unit_E nodes in a_p->u_list if needed and      */
/*   leave b_p->u_list intact                                   */
/*   a_p->u_scale is multiplied by pow(b_p->u_scale,exp_scale)  */
/*  ----------------------------------------------------------- */
void
u_copy_unit(Unit_t *a_p, Unit_t *b_p, double exp_scale) 
{
  Unit_E  *oe_p, *ne_p, *last_p;
  int      ii;
  double   scale;
  
  if( a_p->u_count > 0 ) {
    for(last_p = a_p->u_list; last_p->ue_nextp; last_p = last_p->ue_nextp) {  }
  } else {
    a_p->u_list = last_p = NULL;
  }
  if( b_p->u_count > 0 ) {
    oe_p = b_p->u_list;
    for(ii=0;ii<b_p->u_count;ii++) {
      ne_p = (Unit_E *) capa_malloc(1, sizeof(Unit_E)); /* *** */
      ne_p->ue_scale = oe_p->ue_scale;
      ne_p->ue_exp   = oe_p->ue_exp * exp_scale;
      ne_p->ue_index = oe_p->ue_index;
      strcpy(ne_p->ue_symbol, oe_p->ue_symbol);
      oe_p = oe_p->ue_nextp;
      if( last_p == NULL ) {
        a_p->u_list = ne_p;
      } else {
        last_p->ue_nextp = ne_p;
      }
      last_p = ne_p;
      a_p->u_count++;
    }
    scale = pow(b_p->u_scale, exp_scale);
    a_p->u_scale = a_p->u_scale * scale;
    /* printf("Found scale=%g=%g\n",a_p->u_scale,b_p->u_scale); */
  } else {  
    if( b_p->u_type == U_BASE ) { 
      /* *b_p is a base unit, so create a one element unit */
      ne_p = (Unit_E *) capa_malloc(1, sizeof(Unit_E));   /* *** */
      ne_p->ue_scale = b_p->u_scale;
      ne_p->ue_exp   = exp_scale;
      ne_p->ue_index = b_p->u_index;
      strcpy(ne_p->ue_symbol, b_p->u_symbol);
      if( last_p == NULL ) {
        a_p->u_list = ne_p;
      } else {
        last_p->ue_nextp = ne_p;
      }
      last_p = ne_p;
      a_p->u_count++;
    } else if( b_p->u_type == U_DERIVED) {
      /* derived units but without any units elements (scalar) */
      /*a_p->u_count++;*/
      scale = pow(b_p->u_scale, exp_scale);
      a_p->u_scale = a_p->u_scale * scale;
    } else if( b_p->u_type == U_CONSTANT ) {
      scale = pow(b_p->u_scale, exp_scale);
      a_p->u_scale = a_p->u_scale * scale;
    } else {
      printf("This node has no u_e list and Type unknown\n");
    }
  }
}
int
u_pm_op(Unit_t *a_p, Unit_t *b_p, int op)
{
  int    result=0;
  
  if( a_p->u_count > 0 || b_p->u_count > 0 ) {
     printf(" cannot add or sub units at this moment\n");
     return  result;
  }
  if( op == U_OP_PLUS ) {
    a_p->u_scale = a_p->u_scale + b_p->u_scale;
  } else {
    a_p->u_scale = a_p->u_scale - b_p->u_scale;
  }
  return 1;
}

int
u_parsepower(char *unit_str)
{
  int   exp, ii;
  char  *ch_p, exp_str[16];

  ch_p = unit_str;
  while( isspace(*ch_p) ) { ch_p++; }
  ii=0;
  while( isdigit(*ch_p) ) {
    ch_p++;
  }
  while( isspace(*ch_p) ) { ch_p++; }
  if( *ch_p == '^' ) {
    ch_p++;
  }
  while( isspace(*ch_p) ) { ch_p++; }
  if( *ch_p == '{' ) {
    ch_p++;
  }
  while( isspace(*ch_p) ) { ch_p++; }
  ii=0;
  while( isdigit(*ch_p) || *ch_p == '-' || *ch_p == '+' ) {
    exp_str[ii++] = *ch_p;
    ch_p++;
  }
  exp_str[ii]=0;
  sscanf(exp_str,"%d", &exp);
  return (exp);
}

/* ------------------------------------------- */
/* scan a number of the form indicated below from the input buffer */
/* 1.234^{2.3} */
/*  1e */
double
s_scan_number(char *buf, int idx, int *r_idx)
{
  double   num; 
  float    exp; 
  double   result;
  int      ii=0;
  char     num_str[QUARTER_K];
  
  num_str[ii]=0;
  
  if( buf[idx] == '-' ) {
    num_str[ii++] = '-';
    idx++;
  }
  while( isdigit(buf[idx]) || buf[idx] == '.' ) { 
      num_str[ii++] = buf[idx];
      idx++;
  }
  if( buf[idx] == 'E' || buf[idx] == 'e' ) {
    if( buf[idx+1] == '-' || isdigit(buf[idx+1]) ) {
      num_str[ii++] = buf[idx++];
      num_str[ii++] = buf[idx++];
      while( isdigit(buf[idx]) ) {
        num_str[ii++] = buf[idx];
        idx++;
      }
    }
  }
  num_str[ii] = 0; /* terminate the str */
  sscanf(num_str,"%lg", &num);
  /* printf("Scan number %s got %g\n",num_str, num); fflush(stdout); */
  result = num;
  if( buf[idx] == '^' ) {
    idx++;
    while( isspace(buf[idx]) ) { idx++; }
    if( buf[idx] == '{' ) {  /* need to scan for a matching right bracket */
        idx++;
    }
    while( isspace(buf[idx]) ) { idx++; }
    num_str[0]=0;
    if( isdigit(buf[idx]) || buf[idx] == '+' || buf[idx] == '-' )  {
       ii=0;
       while( isdigit(buf[idx]) || buf[idx] == '.' || buf[idx] == '+' || buf[idx] == '-' ) {
         num_str[ii++] = buf[idx];
         idx++;
       }
       num_str[ii]=0;
    }
    while( isspace(buf[idx]) ) { idx++; }
    if( buf[idx] == '}' ) {
      idx++;
    }
    sscanf(num_str,"%f", &exp);
    /* printf("Scan exp number %s got %g\n",num_str, exp); fflush(stdout); */
    
    result = pow(num, (double)exp);
    /* printf("{%d^%d}=%g\n",num, exp,result); */
  }
  *r_idx = idx;
  return (result);
}


double
s_scan_symbol(char *buf,char *symb_p,int idx, int *r_idx)
{
  char     num_str[QUARTER_K];
  int      ii=0;
  double   r_exp=1.0;
  
  symb_p[0]=0;
  while( isalnum(buf[idx]) || buf[idx] == '_' ) {
    symb_p[ii++] = buf[idx];
    idx++;
  }
  symb_p[ii]=0;
  
  if( buf[idx] == '^' ) {  /* look for either left bracket or a number */
    idx++;
    while( isspace(buf[idx]) ) { idx++; } 
    if( buf[idx] == '{' ) {  /* need to scan for a matching right bracket */
      idx++;
    }
    while( isspace(buf[idx]) ) { idx++; }
    if( isdigit(buf[idx]) || buf[idx] == '.' || buf[idx] == '+' || buf[idx] == '-'  )  {
      ii=0; num_str[ii] = 0;
      while( isdigit(buf[idx]) || buf[idx] == '.' || buf[idx] == '+' || buf[idx] == '-' ) {
        num_str[ii++] = buf[idx];
        idx++;
      }
      num_str[ii]=0;
    }
    while( isspace(buf[idx]) ) { idx++; }
    if( buf[idx] == '}' ) {
      idx++;
    }
    sscanf(num_str,"%lg", &r_exp);  /* power could be of type float */
    /* printf("[scan symb with power %s ^ %lg] ",symb_p, r_exp); fflush(stdout);  */
  }
  *r_idx = idx;
  return (r_exp);
}

/*  return: err_code    0    parsed ok */
/*                      1    symbol is of length 1, not found in the tree */
/*                      2    symbol not found in the tree  */
/*                      3    symbol parsed as prefix symb, but symb not found */
/*                      4    symbol length is 0 or negative */
int
s_process_symb(char *symb_str,Unit_t  *cu_p,double exp)
{
  int      len;
  Unit_t  *au_p;
  int      c_result;
  int      ii;
  char     tmp_str[ANSWER_STRING_LENG];
  int      err_code = 0;
  double   d_exp;
  
  len = strlen(symb_str);
  if( len > 0 ) {
    au_p = u_find_symb(symb_str, UnitTree_p, &c_result);
    if( c_result == 1 ) {  /* if found, copy the definition over */
      u_copy_unit(cu_p, au_p, exp);
    } else {
      if( len > 1 ) {
        if( PrefixTbl[ (int)symb_str[0] ] != 0 ) {  /* prefix is defined */
          for(ii=1;ii<len;ii++) {
             tmp_str[ii-1] = symb_str[ii];
          }
          tmp_str[len-1]=0;
          au_p = u_find_symb(tmp_str, UnitTree_p, &c_result);
          if( c_result == 1 ) {
              /* printf("[%s] ", tmp_str); */
            u_copy_unit(cu_p, au_p, exp);
            d_exp = (double)PrefixTbl[ (int)symb_str[0] ] * exp;
            cu_p->u_scale = cu_p->u_scale * pow((double)10.0,d_exp);
          } else { /* unit *tmp_str not found */
            /*printf("The unit: %s, not defined\n",tmp_str);*/
            err_code = 3;
          }
        } else {
          /*printf("<<%s>>", symb_str);*/
          err_code = 2;
        }
      } else {/* len == 1 */
	/*printf("The unit: %s, not defined\n",symb_str);*/
        err_code = 1;
      }
    }
  } else {
    err_code = 4;
  }
  return (err_code);
}

Unit_t *
u_parse_unit(char *unit_str)
{
  char      *ch;
  char       symb_str[QUARTER_K];
  int        idx;
  double     exp_sign;
  int        s_result;
  int        not_done;
  double     s_number,  offset;
  double     tmp_scale, symb_exp, exp;
  Unit_t    *cu_p;

  gUnitError=0;
  ch   = unit_str;
  cu_p = (Unit_t *) capa_malloc(1, sizeof(Unit_t)); /* *** */
  cu_p->u_scale = 1.0;
  idx = 0;  not_done = 1;
  exp_sign = 1.0; exp = 1;
  symb_str[0] = 0;

  while( isspace(*ch) ) { ch++; }    /* trim leading white spaces */
  /* fprintf(stdout,"PARSE |%s|\n", unit_str); */
  while( not_done ) {
    if( isdigit(ch[idx]) || ch[idx] == '-' ) {  /* rule 1: number */
       s_number = s_scan_number(ch,idx,&idx);
       
       tmp_scale = pow(s_number,exp_sign);
       /* printf("S=%g,Power(%g,%d)=%g\n", 
          cu_p->u_scale, s_number,exp_sign, tmp_scale);
       */
       cu_p->u_scale = cu_p->u_scale * tmp_scale;
       
       /* printf("[Scale %g=%g^%g] ",tmp_scale,s_number,exp_sign); */
       while( isspace(ch[idx]) ) { idx++; }
    } else {
      if( isalpha(ch[idx]) ) { /* rule 2: unit_symbol ^ exp */
	symb_str[0] = 0;
	symb_exp = s_scan_symbol(ch,symb_str,idx,&idx);
	exp = (double)exp_sign * symb_exp;
	/* printf("[scanned %s ^ (%g * %g)] ", symb_str,symb_exp,exp_sign); fflush(stdout); */
	s_result = s_process_symb(symb_str,cu_p,exp);
	if( s_result > 0 ) {
	  /* printf("Error processing symbol [%s]\n", symb_str); */
	  gUnitError = 1;
	}
	while( isspace(ch[idx]) ) { idx++; }
      } else {
	if( ch[idx] == '*' || ch[idx] == '/' ) {
	  if( ch[idx] == '/' ) { /* printf("[/] "); */ exp_sign = -1.0; }
	  idx++;
	  while( isspace(ch[idx]) ) { idx++; }
	} else {
	  if( ch[idx] == '+' || ch[idx] == '-' ) {
	    idx++;
	    while( isspace(ch[idx]) ) { idx++; }
	    offset = s_scan_number(ch,idx,&idx);
	    /* printf("[Offset %g] ",offset); */
	  } else {
	    if( ch[idx] == 0 ) {  /* end of input string */
	      not_done = 0;
	      /* printf("\n"); */
	    } else {
	      /* garbage in unit string */
	      gUnitError = 1;
	      not_done=0;
	    }
	  }
	}
      }
    }
  }
  simplify_unit(cu_p);
  return (cu_p);

}

void
u_getunit(FILE *f)
{
  register int  unit_type;
  register int  c;
  int      power, result;
  char   *name_p, *symbol_p, *comment_p, *unit_p;
  
  BaseUnitcnt = 0;
  free_utree(UnitTree_p);
  UnitTree_p = NULL;
  c_moveto_unit(f);  /* move the file position to << */
  do {
    c_ignorewhite(f);
    c = getc(f); ungetc(c,f);
    if( c == '<' ) {
      unit_type = c_gettype(f);
    }
    if( c != EOF ) {
      switch(unit_type) {
        case U_BASE:
               name_p    = c_getword(f);    symbol_p = c_getword(f); 
               comment_p = c_getcomment(f);
               /*
               printf("B Unit: N=%s,S=%s,C=%s\n",name_p,symbol_p,comment_p);
               */
               result = u_insert_baseunit(name_p,symbol_p,comment_p);
               if( result == 1 ) {
                 printf("The entry %s is duplicated\n",symbol_p);
               }
               free(name_p); free(symbol_p); free(comment_p);
               break;
        case U_DERIVED:
               name_p    = c_getword(f);    symbol_p = c_getword(f);
               unit_p    = c_getstring(f);  comment_p = c_getcomment(f);
               /*
               printf("D Unit: N=%s,S=%s,C=%s,U=%s\n",
                       name_p,symbol_p,comment_p,unit_p);
               */
               result = u_insert_derived(name_p,symbol_p,comment_p,unit_p);
               if( result == 1 ) {
                 printf("The entry %s is duplicated\n",symbol_p);
               }
               /* preorder_utree(UnitTree_p); */ 
               free(name_p); free(symbol_p); free(comment_p); free(unit_p);
               break;
        case U_PREFIX:
               name_p    = c_getword(f);    symbol_p = c_getword(f);
               unit_p    = c_getstring(f);
               /*
               printf("Prefix: N=%s,S=%s,U=%s\n",
                       name_p,symbol_p,unit_p);
               */
               power = u_parsepower(unit_p);
               PrefixTbl[ (int)(*symbol_p) ] = power;
               /* printf("    P[%c]=%d\n",*symbol_p,power);  */
               free(name_p); free(symbol_p); free(unit_p);
               break;
        case U_CONSTANT:
               symbol_p = c_getword(f);  unit_p    = c_getstring(f);
               comment_p = c_getcomment(f);
               /*
               printf("Const.: S=%s,C=%s,U=%s\n",
                       symbol_p,comment_p,unit_p);
               */
               break;
        case U_UNKNOWN:
               /* printf("Unknown\n"); */
               break;
      }
    }
  } while ( c != EOF );

}

/* ----------------------------------------------------------------- */
/* comparing unit symbol names should be case sensitive */
int
comp_unit_symb(a, b) char *a; char *b;
{
  return strncmp(a,b,SYMBOL_MAXLEN);
}


Unit_t *
u_splay (char *name, Unit_t *t) 
{
  Unit_t     N;
  Unit_t    *l, *r, *y;

  if (t == NULL)  return t;
  N.u_left  = (Unit_t *)NULL;
  N.u_right = (Unit_t *)NULL;
  l = r = &N;

  for (;;) {
    if ( comp_unit_symb(name,t->u_symbol) < 0 ) {
      if (t->u_left == NULL)  break;
      if ( comp_unit_symb(name, (t->u_left)->u_symbol ) < 0 ) {
        y = t->u_left; t->u_left = y->u_right; y->u_right = t; t = y;
        if (t->u_left == NULL) break;
      }
      r->u_left = t; r = t; t = t->u_left;
    } else if ( comp_unit_symb(name,t->u_symbol) > 0 ) {
        if (t->u_right == NULL) break;
        if ( comp_unit_symb(name, (t->u_right)->u_symbol ) > 0 ) {
          y = t->u_right; t->u_right = y->u_left; y->u_left = t; t = y;
          if (t->u_right == NULL) break;
        }
        l->u_right = t; l = t; t = t->u_right;
    } else {
      break;
    }
  }
  l->u_right = t->u_left; r->u_left = t->u_right; t->u_left = N.u_right;
  t->u_right = N.u_left;
  return t;
}



/* returns: 0  correctly inserted */
/*          -1 error */
/*          1  duplicate entry    */

int
u_insert_baseunit(n_p,s_p,c_p) char  *n_p, *s_p, *c_p;
{
  Unit_t   *new_p, *t;
  int       len;
 
  new_p = (Unit_t *) capa_malloc(1, sizeof(Unit_t)); /* *** */
  if (new_p == NULL) {
      printf("Ran out of space\n");
      return(-1);
  }
  strcpy(new_p->u_symbol, s_p);
  strcpy(new_p->u_name, n_p);
  len = strlen(c_p);
  new_p->u_comment = (char *) capa_malloc((len+1), sizeof(char)); /* *** */
  strcpy(new_p->u_comment,c_p);
  BaseUnitcnt++;
  new_p->u_index  = BaseUnitcnt;
  new_p->u_type   = U_BASE;
  new_p->u_scale  = 1.0;
  new_p->u_offset = 0.0;
  new_p->u_count  = 0;
  new_p->u_list   = NULL;
 
  if (UnitTree_p == NULL) {  /* a new unit tree */
      UnitTree_p = new_p;
      return (0);
  }
  t = u_splay(s_p, UnitTree_p);
  if ( comp_unit_symb(s_p,t->u_symbol) < 0 ) {
        new_p->u_left = t->u_left; new_p->u_right = t;
        t->u_left = NULL;
        /* Splay_cnt++;  */
        UnitTree_p = new_p;
        return (0);
  } else if ( comp_unit_symb(s_p,t->u_symbol) > 0 ) {
        new_p->u_right = t->u_right; new_p->u_left = t;
        t->u_right = NULL;
        /* Splay_cnt++; */
        UnitTree_p = new_p;
        return (0);
  } else {    /* name and t->u_symbol is the same, which means found it */
        capa_mfree( (char *)new_p );
        UnitTree_p = t;
        return (1);
  }
}


int
u_insert_derived(n_p,s_p,c_p,u_p)char  *n_p, *s_p, *c_p, *u_p;
{
  Unit_t  *new_p, *t;
  int      c_result, len;
  
  /* inorder_utree(UnitTree_p); */
  t = u_splay(s_p, UnitTree_p);
  UnitTree_p = t;
  c_result = comp_unit_symb(s_p,t->u_symbol);
  if ( c_result == 0 ) {
    UnitTree_p = t;
    return (1);
  }
  
  /* prepare a new Unit_t */
  new_p = u_parse_unit(u_p);
  strcpy(new_p->u_symbol,s_p);
  strcpy(new_p->u_name, n_p);
  new_p->u_type = U_DERIVED;
  len = strlen(c_p);
  new_p->u_comment = (char *) capa_malloc((len+1), sizeof(char)); /* *** */
  strcpy(new_p->u_comment,c_p);
  
  simplify_unit(new_p);
#ifdef UNIT_DBUG
  printf("Derived Unit:%s\n",new_p->u_name);
  print_unit_t(new_p); 
#endif
  if (c_result < 0 ) {
    new_p->u_left = t->u_left; new_p->u_right = t;
    t->u_left = NULL;
  } else {  /* c_result > 0 */
    new_p->u_right = t->u_right; new_p->u_left = t;
    t->u_right = NULL;
  }
  UnitTree_p = new_p;
  
  return (0);
  
}

void
freelist_unit_e(Unit_E *ue_p) 
{
  Unit_E  *curr_p, *next_p;
  
  if( ue_p != NULL ) {
    next_p = ue_p->ue_nextp;
    curr_p = ue_p;
    if( next_p == NULL ) {
      capa_mfree((char *)curr_p);
    } else {
      for( curr_p = ue_p; next_p; curr_p = next_p, next_p = next_p->ue_nextp) {
        capa_mfree((char *)curr_p);
      }
      capa_mfree((char *)curr_p);
    }
  }
}
void
simplify_unit(u_p) Unit_t *u_p;
{
  Unit_E   *eu_p, *prev_p;
  int       ii, idx;
  
  /* walk through u_list and replace those u_index = -1 with */
  /* a linked list of basic unit. */
  /* u_msort_main() the whole u_list */
  /* combine those units with same u_index */
  for(ii=0;ii<BaseUnitcnt;ii++) {
    CScale[ii] = 0.0;
    CExp[ii] = 0.0;
  }
  /*
  printf("Before Simplify:: \n");
  print_unit_t(u_p);
  */
  if( u_p->u_count > 0 ) {
    
    for(eu_p=u_p->u_list; eu_p; eu_p = eu_p->ue_nextp) {
      idx = eu_p->ue_index;
      if( CScale[idx] == 0.0 ) {
        CScale[idx] = 1.0;
        strcpy(CSymb[idx],eu_p->ue_symbol);
      }
      CScale[idx] = CScale[idx] * eu_p->ue_scale;
      CExp[idx] = CExp[idx] + eu_p->ue_exp;
    }
    /* debugging 
    for(ii=0;ii<BaseUnitcnt;ii++) {
      if( CScale[ii] != 0.0 ) {
        printf("(%d)%s,S=%g,E=%g\n",ii,CSymb[ii],CScale[ii], CExp[ii]);
      }
      if( CExp[ii] == 0.0 ) {
        printf("(%d)%s,S=%g,Exp=%g\n",ii,CSymb[ii],CScale[ii], CExp[ii]);
      }
    }
    */
    freelist_unit_e(u_p->u_list);
    prev_p = u_p->u_list = NULL;
    u_p->u_count = 0;
    for(ii=0;ii<BaseUnitcnt;ii++) {
      if( CScale[ii] != 0.0 && CExp[ii] != 0) {
        eu_p = (Unit_E *)capa_malloc(1,sizeof(Unit_E)); /* ***************** */
        eu_p->ue_scale = 1.0;
        eu_p->ue_exp = CExp[ii];
        eu_p->ue_index = ii;
        strcpy(eu_p->ue_symbol,CSymb[ii]);
        if( prev_p == NULL) {
          u_p->u_list = prev_p = eu_p;
        } else {
          prev_p->ue_nextp = eu_p;
          prev_p = eu_p;
        }
        u_p->u_count++;
      }
    }
  }
  /* 
  printf("After Simplify:: \n");
  print_unit_t(u_p);
  */
}

/* before comparing two units, make sure they are of  basic form */
/* compares if two units are equal */
/* equality returns 1 */

int  is_units_equal(Unit_t *u1_p, Unit_t *u2_p)
{
  int      result=1;
  Unit_E  *a_p, *b_p;
  
  if( (u1_p->u_count == u2_p->u_count) && 
      (u1_p->u_scale == u2_p->u_scale) ) {
    for(a_p=u1_p->u_list, b_p=u2_p->u_list;
        a_p; a_p=a_p->ue_nextp, b_p=b_p->ue_nextp) {
      if(a_p->ue_index != b_p->ue_index ||
         a_p->ue_scale != b_p->ue_scale ||
         a_p->ue_exp   != b_p->ue_exp ) {
        result=0;
        break;
      }
    }
  } else {
    result=0;
  }
  return (result);
}
/*     input : both are the simplest units */
/*     result: 0.0 means they are not of euquvalent units */
/*             the ratio of u1 / u2   */
double  units_ratio(Unit_t *u1_p, Unit_t *u2_p)
{
  double   ratio=1.0;
  Unit_E  *a_p, *b_p;
  
  if( (u1_p->u_count == u2_p->u_count) ) {
    for(a_p=u1_p->u_list, b_p=u2_p->u_list;
        a_p; a_p=a_p->ue_nextp, b_p=b_p->ue_nextp) {
      if(a_p->ue_index != b_p->ue_index ||
         a_p->ue_scale != b_p->ue_scale ||
         a_p->ue_exp   != b_p->ue_exp ) {
        ratio=0.0;
        break;
      }
    }
  } else {
    ratio=0.0;
  }
  if( (ratio != 0.0) && (u2_p->u_scale != 0.0 )  ) {
    ratio = u1_p->u_scale / u2_p->u_scale;
  }
  return (ratio);
}

/* ------------- The Grammar of Units Parser --------------------

  scan_unit_expr()  -->  scan_basic_block()
                    -->  scan_basic_block() '+' scan_basic_block() 
                    -->  scan_basic_block() '-' scan_basic_block()
 
  scan_num_expr()   -->  scan_num_block()
                    -->  scan_num_block() '+' scan_num_block()
                    -->  scan_num_block() '-' scan_num_block()
                    
  scan_basic_block()-->  scan_basic_term()
                    -->  scan_basic_term()  '*' scan_basic_term()
                    -->  scan_basic_term()  ' ' scan_basic_term()
                    -->  scan_basic_term()  '/' scan_basic_term()

  scan_num_block()  -->  scan_num_term()
                    -->  scan_num_term()  '*' scan_num_term()
                    -->  scan_num_term()  ' ' scan_num_term()
                    -->  scan_num_term()  '/' scan_num_term()
  
  
  scan_basic_term() -->  scan_unit_item()          
                    -->  scan_num_item()
                    -->  '(' scan_basic_block() ')'
                    -->  '{' scan_basic_block() '}'

  scan_num_term()   -->  scan_num_item()<sp>*
                    --> '-' scan_num_item()<sp>*
                    --> '(' scan_num_expr() ')'
                    --> '{' scan_num_expr() '}'

  scan_unit_item()  -->  UNIT<sp>*
                    -->  UNIT<sp>*  '^' <sp>* scan_num_term()
                    
  scan_num_item()   -->  FLOAT<sp>*
                    -->  FLOAT<sp>* '^' <sp>* scan_num_term()
  
  scan_FLOAT()      -->  [0-9]+([eE][+-]?[0-9]+)*
  
  p_new_unit()      -->  [a-Z]+[a-Z0-9_]*
  
  -----------------------------------------
  U.expr  := B.block
           | B.block '+' B.block
           | B.block '-' B.block
           
  N.expr  := N.block 
           | N.block '+' N.block
           | N.block '-' N.block
 
 To allow for operations like (J/N)^2 or {N/m}^2 (N/J)^3 
 
 
  B.block := B.term
           | B.term ' ' B.term
           | B.term '*' B.term
           | B.term '/' B.term
           
  N.block := N.term 
           | N.term ' ' N.term
           | N.term '*' N.term
           | N.term '/' N.term
           
  B.term  := U.item
           | N.item
           | '(' B.block ')'
           | '{' B.block '}'
           
           | '(' B.block ')' ^ N.term
           | '{' B.block '}' ^ N.term
           
  N.term  := N.item
           | '-' N.item
           | '(' N.expr ')'
           | '{' N.expr '}'
           
  U.item  := UNIT
           | UNIT '^' N.term
           
  N.item  := FLOAT
           | FLOAT '^' N.term
           
  UNIT    := [a-Z]+[a-Z0-9_]*
  
  FLOAT   := [0-9]+([eE][+-]?[0-9]+)*
  
 ------------------------------------------------------------------- */
 
Unit_t *
p_new_op(Unit_t *left_p, int op, Unit_t *right_p)
{
  Unit_t  *new_p;
  
  new_p = (Unit_t *) capa_malloc(1, sizeof(Unit_t));
  if (new_p == NULL) {
      printf("Ran out of space\n");
      return(NULL);
  }
  new_p->u_left   = left_p;
  new_p->u_right  = right_p;
  new_p->u_scale  = 0.0;
  new_p->u_type   = op;
  new_p->u_offset = 0.0;
  new_p->u_count  = 0;
  new_p->u_list   = NULL;
  
  return (new_p);
}

Unit_t *
p_new_num(Unit_t *left_p, double num, Unit_t *right_p)
{
  Unit_t  *new_p;
  
  new_p = (Unit_t *) capa_malloc(1, sizeof(Unit_t));
  if (new_p == NULL) {
      printf("Ran out of space\n");
      return(NULL);
  }
  
  new_p->u_left   = left_p;
  new_p->u_right  = right_p;
  new_p->u_scale  = num;
  new_p->u_type   = U_CONSTANT;
  new_p->u_offset = 0.0;
  new_p->u_count  = 0;
  new_p->u_list   = NULL;
  
  return (new_p);
}

Unit_t *
p_new_unit(Unit_t *left_p, Unit_t *right_p)
{
  char     symb_str[ANSWER_STRING_LENG];
  int      ii=0;
  int      len;
  Unit_t  *au_p, *cu_p;
  int      c_result;
  char     tmp_str[ANSWER_STRING_LENG];
  int      err_code = 0;
  double   d_exp;
  
  symb_str[ii]=0;
  while( isspace(Sbuf[Sidx]) ) { Sidx++; }
  while( isalnum(Sbuf[Sidx]) || Sbuf[Sidx] == '_' ) {
    symb_str[ii++] = Sbuf[Sidx];
    Sidx++;
  }
  symb_str[ii]=0;
  /* printf("<U %s>", symb_str); */
  cu_p = (Unit_t *) capa_malloc(1, sizeof(Unit_t));
  strcpy(cu_p->u_symbol,symb_str);
  cu_p->u_left   = left_p;
  cu_p->u_right  = right_p;
  cu_p->u_scale  = 1.0;
  cu_p->u_type   = U_DERIVED;
  cu_p->u_offset = 0.0;
  cu_p->u_count  = 0;
  cu_p->u_list   = NULL;
  
  len = strlen(symb_str);
  if( len > 0 ) {
    au_p = u_find_symb(symb_str, UnitTree_p, &c_result);
    if( c_result == 1 ) {  /* if found, copy the definition over */
      u_copy_unit(cu_p, au_p, 1);
    } else {
      if( len > 1 ) {
        if( PrefixTbl[ (int)symb_str[0] ] != 0 ) {  /* prefix is defined */
          for(ii=1;ii<len;ii++) {
             tmp_str[ii-1] = symb_str[ii];
          }
          tmp_str[len-1]=0;
          au_p = u_find_symb(tmp_str, UnitTree_p, &c_result);
          if( c_result == 1 ) {
              /* printf("[%s] ", tmp_str); */
            u_copy_unit(cu_p, au_p, 1);
            d_exp = (double)PrefixTbl[ (int)symb_str[0] ];
            cu_p->u_scale = cu_p->u_scale * pow((double)10.0,d_exp);
          } else { /* unit *tmp_str not found */
            /* printf(" not found\n"); */
            err_code = 3;
	    cu_p->u_type   = U_UNKNOWN;
          }
        } else { /* symb_str is not in <prefix><units> form */
          /* printf("<<%s>>", symb_str); */
          err_code = 2;
	  cu_p->u_type   = U_UNKNOWN;
        }
      } else {/* len == 1 */
        /* printf(" not found in symbol tree \n"); */
        err_code = 1;
	cu_p->u_type   = U_UNKNOWN;
      }
    }
  } else { /* why would we have a length less than zero symb_str ? */
    err_code = 4;
  }
  
  return (cu_p);
}

int  s_peeknext_op()
{
  char  *ch;
  int    sp=0;
  
  ch = (char *)&Sbuf[Sidx];
  while( isspace(*ch) ) { ch++; sp=1; }
  if( (*ch == '*')  || (*ch == '/') || (*ch == '+')  || (*ch == '-') || (*ch == '^')) {
    return (*ch);
  }
  /* what if space is the last thing on the line?*/
  if( sp && (*ch != '\0')) return '*';
  return (*ch);
}

int  s_getnext_op()
{
  char  *ch;
  int    inc = 0, sp=0;
  
  
  /* printf("\n((op"); print_remains(); printf("\n");  */
  ch = (char *)&Sbuf[Sidx];
  while( isspace(*ch) ) { ch++; inc++; sp=1; }
  Sidx = Sidx + inc;
  if( (*ch == '*')  || (*ch == '/') || (*ch == '+')  || (*ch == '-') || (*ch == '^') ) {
    Sidx++;
    /* print_remains();  printf(" op))"); printf("\n"); */
    return (*ch);
  }
  /* print_remains();  printf(" op))"); printf("\n"); */
  /* what if space is the last thing on the line?*/
  if( sp  && (*ch != '\0')) return '*';
  return (*ch);
}

int
s_getnext()
{
  char  ch;
  
  ch = Sbuf[Sidx];
  Sidx++;
  return (ch);
}

int
s_peeknext()
{
  char  ch;
  
  ch = Sbuf[Sidx];
  return (ch);
}

int
s_peeknextNW()  /* peek into the next non-whitespaces character */
{
  char  *ch;

  ch = (char *)&Sbuf[Sidx];
  while( isspace(*ch) ) { ch++; }
  return (*ch);
}

int
s_getnextNW()  /* get the next non-whitespaces character */
{
  char  *ch;

  ch = (char *)&Sbuf[Sidx]; Sidx++;
  while( isspace(*ch) ) { ch++; Sidx++; }
  return (*ch);
}
/* peek into the next non-whitespaces character 
   which should be either a multiply or division */
int
s_peekMDWS()  
{
  char  *ch;
  int    sp=0;
  
  ch = (char *)&Sbuf[Sidx];
  while( isspace(*ch) ) { ch++; sp=1;}
  if( (*ch == '*')  || (*ch == '/') ) {
    return (*ch);
  }
  if( sp ) return ' ';
  ch = (char *)&Sbuf[Sidx];
  while( isspace(*ch) ) { ch++; }
  return (*ch);
}

int
s_getnextMDWS()
{
  char  *ch;
  int    inc=0, sp=0;
  
  ch = (char *)&Sbuf[Sidx]; Sidx++;
  while( isspace(*ch) ) { ch++; inc++; sp=1; }
  Sidx += inc;
  if( (*ch == '*')  || (*ch == '/') ) {
    return (*ch);
  }
  if( sp ) return ' ';
  return (*ch);
}

double
scan_FLOAT()
{
  double   num; 
  int      ii=0, len;
  char     num_str[QUARTER_K];
  
  num_str[ii]=0;
  while( isspace(Sbuf[Sidx]) ) { Sidx++; }
  if( Sbuf[Sidx] == '-' ) {
    num_str[ii++] = Sbuf[Sidx++];
  }
  while( isdigit(Sbuf[Sidx]) || Sbuf[Sidx] == '.' ) {
      num_str[ii++] = Sbuf[Sidx++];
  }
  if( Sbuf[Sidx] == 'E' || Sbuf[Sidx] == 'e' ) {
    if( Sbuf[Sidx+1] == '-' || isdigit(Sbuf[Sidx+1]) ) {
      num_str[ii++] = Sbuf[Sidx++];
      num_str[ii++] = Sbuf[Sidx++];
      while( isdigit(Sbuf[Sidx]) ) {
        num_str[ii++] = Sbuf[Sidx++];
      }
    }
  }
  num_str[ii] = 0; /* terminate the str */
  len = strlen(num_str);
  if(len > 0 ) {
    sscanf(num_str,"%lg", &num);
    /* printf("<N %s %g>",num_str,num); fflush(stdout);  print_remains(); */
  } else {
    num = 1.0;
  }
  return (num);
}
/* -----------------------------------------------
  N.item  := FLOAT
           | FLOAT '^' N.term
   ----------------------------------------------- */
Unit_t  *
scan_num_item()
{
  Unit_t  *node_p, *exp_p;
  double   num_const;
  char     ch;
  
  num_const = scan_FLOAT();
  node_p = p_new_num(NULL, num_const, NULL);
  ch = s_peeknext_op();
  if( ch == '^' ) {
    ch = s_getnext_op();
    
    exp_p = scan_num_term();
    num_const = node_p->u_scale;
    if( node_p->u_scale > 0.0 ) {
      num_const = pow(node_p->u_scale,exp_p->u_scale);
    }
    node_p->u_scale = num_const;
    capa_mfree((char *)exp_p);
  }
  return node_p;
}

/* -----------------------------------------------
  U.item  := UNIT
           | UNIT '^' N.term
   ----------------------------------------------- */
   
Unit_t *
scan_unit_item()
{
  Unit_t   *node_p, *exp_p;
  char      ch;
  double   num_const;
  Unit_E   *oe_p;
  
  node_p = p_new_unit(NULL,NULL);
  ch = s_peeknext_op();
  if( ch == '^' ) {
    ch = s_getnext_op();
    exp_p = scan_num_term();
    num_const = exp_p->u_scale;
    if( node_p->u_count > 0 ) {
      oe_p = node_p->u_list;
      for(oe_p = node_p->u_list; oe_p; oe_p = oe_p->ue_nextp ) {
        oe_p->ue_exp   = oe_p->ue_exp * num_const;
      }
    }
    num_const = node_p->u_scale;
    if( node_p->u_scale > 0.0 ) {
      num_const = pow(node_p->u_scale,exp_p->u_scale);
    }
    node_p->u_scale = num_const;
    capa_mfree((char *)exp_p);
  }
  return node_p;
}

void distribute_exp(Unit_t* node_p,Unit_t* exp_p) 
{
  Unit_E* oe_p;
  double num_const;
  num_const = exp_p->u_scale;  /* should we check if num_const too large or small ? */
  if( node_p->u_count > 0 ) {
    oe_p = node_p->u_list;
    for(oe_p = node_p->u_list; oe_p; oe_p = oe_p->ue_nextp ) {
      oe_p->ue_exp   = oe_p->ue_exp * num_const;
    }
  }
  num_const = node_p->u_scale;
  if( node_p->u_scale > 0.0 ) {  /* what if u_scale <= 0.0 ? */
    num_const = pow(node_p->u_scale,exp_p->u_scale);
  }
  node_p->u_scale = num_const;
  if (node_p->u_left) distribute_exp(node_p->u_left,exp_p);
  if (node_p->u_right) distribute_exp(node_p->u_right,exp_p);
}

/* ---------------------------------------------------------------
   B.term  := U.item
           | N.item
           | '(' B.block ')'
           | '{' B.block '}'
           
           | '(' B.block ')' '^' N.term  <== July 6 1998
           | '{' B.block '}' '^' N.term
           
   --------------------------------------------------------------- */
Unit_t *
scan_basic_term()
{
  Unit_t   *node_p, *exp_p;
  int       ch, nch;
  
  ch = s_peeknextNW();
  if( ch == '(' || ch == '{' ) {
    ch = s_getnextNW();  /* get rid of '(' or '{' */
    node_p = scan_basic_block();
    nch = s_peeknextNW();
    if( nch == ')' || nch == '}' ) {  /* should be either ')' or '}' */
      if( ((ch == '(' ) && (nch == ')' )) ||
          ((ch == '{' ) && (nch == '}' )) ) { /* matching left paren with right paren */
          
           
      } else {
        /* printf(" WARN: %c matched by %c\n", ch, nch); */
      }
      nch = s_getnextNW();
      /* ====== Added Jul 6, 1998 ====> */
      ch = s_peeknext_op();
      if( ch == '^' ) {
        ch = s_getnext_op();  /* get rid of '^' char */
        exp_p = scan_num_term();
	distribute_exp(node_p,exp_p);
        capa_mfree((char *)exp_p);
      } 
      /* <== added Jul 6, 1998 == */
    } else {
      /* printf(" WARN: %c is not matched by %c\n", ch, nch); */
    }
  } else if( ch >= '0' && ch <= '9' ) {
    node_p = scan_num_item();
  } else { /* assume a unit symbol */
    /* printf("<B.term>"); print_remains(); */
    node_p = scan_unit_item();
    /* print_remains(); */
  }
  return node_p;
}
/* --------------------------------------------------
   N.term  := N.item
           | '-' N.item
           | '(' N.expr ')'
           | '{' N.expr '}'
 -------------------------------------------------- */
Unit_t *
scan_num_term()
{
  Unit_t   *node_p;
  char      ch, nch;

  ch = s_peeknextNW();
  if( ch == '(' || ch == '{' ) {
    ch = s_getnextNW();
    node_p = scan_num_expr();
    nch = s_peeknextNW();
    if( nch == ')' || nch == '}' ) {  /* should be either ')' or '}' */
      if( ((ch == '(' ) && (nch == ')' )) ||
          ((ch == '{' ) && (nch == '}' )) ) { 
        
      } else {
        /* printf(" WARN: %c matched by %c\n", ch, nch); */
      }
      nch = s_getnextNW();
    } else {
      /* printf(" WARN: %c is not matched by %c\n", ch, ch); */
    }
  } else if( ch == '-' ) {
    ch = s_getnextNW();
    node_p = scan_num_item();
    node_p->u_scale = (-1)*node_p->u_scale;
  } else {
    if( isdigit(ch) ) {
       node_p = scan_num_item();
    } else { /* something other than a number */
       /*
          printf(" ERROR: expect a number: ");
          print_remains();
       */
       node_p = p_new_num(NULL, 0.0, NULL); /* make the unknown item */
    }
  }
  return node_p;
}

/* --------------------------------------------------
   B.block := B.term
           | B.term ' ' B.term
           | B.term '*' B.term
           | B.term '/' B.term
   -------------------------------------------------- */
Unit_t  *
scan_basic_block()
{
  Unit_t   *node_p;
  char      ch;
  int       op;
  
  /* printf("<B.block>(before B.term)"); print_remains(); */
  node_p = scan_basic_term();
  ch = s_peeknext_op();
  while ( ch == '*' || ch == '/' ) {
    op = ( ch == '/' ? U_OP_DIVIDE : U_OP_TIMES);
    ch = s_getnext_op();
    /* printf("<B.block>(/ *)"); print_remains();  */
    node_p = p_new_op(node_p,op,scan_basic_term());
    ch = s_peeknext_op();
  }
  return node_p;
}
/* --------------------------------------------------
   N.block := N.term 
           | N.term ' ' N.term
           | N.term '*' N.term
           | N.term '/' N.term
   -------------------------------------------------- */
Unit_t  *
scan_num_block()
{
  Unit_t   *node_p, *opand_p;
  char      ch;
  double    result;
  
  node_p = scan_num_term();
  ch = s_peeknext_op();
  while ( ch == '*' || ch == '/' ) {
    s_getnext_op();
    opand_p = scan_num_term();
    if( ch == '*' ) {
      result = node_p->u_scale * opand_p->u_scale;
    } else {
      result = node_p->u_scale / opand_p->u_scale;
    }
    node_p->u_scale = result;
    capa_mfree((char *)opand_p);
    ch = s_peeknext_op();
  }
  return node_p;
}

/* ---------------------------------------
   U.expr  := B.block
           | B.block '+' B.block
           | B.block '-' B.block
   --------------------------------------- */
Unit_t  *
scan_unit_expr()
{
  Unit_t   *node_p;
  char      ch;
  int       op;
  
  /* printf("<U.expr>"); print_remains();  */
  node_p = scan_basic_block();
  ch = s_peeknext_op();
  while ( ch == '+' || ch == '-' ) {
    op = ( ch == '+' ? U_OP_PLUS : U_OP_MINUS);
    ch = s_getnext_op();
    /* printf("<U.expr>(+-)"); print_remains(); */
    node_p = p_new_op(node_p,op,scan_basic_block());
    ch = s_peeknext_op();
  }
  return node_p;
}
/* -----------------------------------------
   N.expr  := N.block 
           | N.block '+' N.block
           | N.block '-' N.block
   ----------------------------------------- */
Unit_t  *
scan_num_expr()
{
  Unit_t   *node_p, *opand_p;
  char      ch;
  double    result;
  
  node_p = scan_num_block();
  ch = s_peeknext_op();
  while ( ch == '+' || ch == '-' ) {
    ch = s_getnext_op();
    opand_p = scan_num_block();
    if( ch == '+' ) {
      result = node_p->u_scale + opand_p->u_scale;
    } else {
      result = node_p->u_scale - opand_p->u_scale;
    }
    node_p->u_scale = result;
    capa_mfree((char *)opand_p);
    ch = s_peeknext_op();
  }
  return node_p;
}

/* ----------------------------------------------------------------------- */
/* <--  This is the major entry point to parse an units expression ------> */
Unit_t  *
parse_unit_expr(char *symb_str)
{
  Unit_t   *root_p;
  int       len;
  
  len = strlen(symb_str);
  strcpy(Sbuf,symb_str);  /* copy it into the global Sbuf */
  Sidx=0;
  root_p = scan_unit_expr();
  if(Sidx < len-1 ) {
    /* printf(" WARN: NOT PARSED:");  print_remains(); */
  }
  return (root_p);

}

void
print_remains()
{
  int       len, ii;
  
  len = strlen(Sbuf);
  printf("[[");
  for(ii=Sidx;ii<len;ii++) {
      printf("%c",Sbuf[ii]);
  }
  printf("]]");
  
}



/* =================================================================== */