# include "opihi.h"

/* the result of a matrix operation must go into a temp variable,
   labeled by "m". if one of the input matrices is already a temp variable, we 
   can continue using it this round 
   */

int VV_binary (StackVar *OUT, StackVar *V1, StackVar *V2, char *op) {

  int i, Nx;
  float *out, *M1, *M2;
  char line[512]; // this is only used to report an error 
  
  Nx = V1[0].vector[0].Nelements;
  if (Nx != V2[0].vector[0].Nelements) {
    return (FALSE);
  }

  if (V1[0].type == 'v') {  /** use V1 as temp buffer **/
    OUT[0].vector = V1[0].vector;
    V1[0].type = 'V'; /* prevent it from being freed below */
  } else {
    if (V2[0].type == 'v') { /** use V2 as temp buffer, but header of V1 **/
      OUT[0].vector = V2[0].vector;
      V2[0].type = 'V'; /* prevent it from being freed below */
    } else {  /* no spare temp buffer */
      OUT[0].vector = InitVector ();
      CopyVector (OUT[0].vector, V1[0].vector);
    }
  }
  OUT[0].type = 'v'; /*** <<--- says this is a temporary matrix ***/

  M1  = V1[0].ptr;
  M2  = V2[0].ptr;
  out = OUT[0].ptr = (float *)OUT[0].vector[0].elements;

// XXX this is not quite there:
// 1) do I drop the .ptr element of StackVar?
// 2) things below are using type from a StackVar to get the vector type -- wrong place
// 3) the temp would need lots of work -- is it worth it?
// 4) some operators force the output to be float -- we cannot use a single bit of code to 
//    define the output type based on which is a float and which is not.  getting ugly.
// ??? how much more work to bring the MV, etc functions into shape ???

# define VV_FUNC(OP) \
  if ((M1->type == OPIHI_FLOAT) && (M2->type == OPIHI_FLOAT)) { \
    // output must be float, can use either as a temp \
    opihi_float *M1 = V1[0].vector[0].elements; \
    opihi_float *M2 = V2[0].vector[0].elements; \
    for (i = 0; i < Nx; i++, out++, M1++, M2++) { \
  	*out = *M1 OP *M2; \
    } \
    break; \
  } \
  if ((M1->type != OPIHI_FLOAT) && (M2->type == OPIHI_FLOAT)) { \
    // output must be float, can use M2 as a temp \
    opihi_float *M1 = V1[0].vector[0].elementsInt; \
    opihi_float *M2 = V2[0].vector[0].elements; \
    for (i = 0; i < Nx; i++, out++, M1++, M2++) { \
  	*out = *M1 OP *M2; \
    } \
    break; \
  } \
  if ((M1->type == OPIHI_FLOAT) && (M2->type != OPIHI_FLOAT)) { \
    // output must be float, can use M1 as a temp \
    opihi_float *M1 = V1[0].vector[0].elements; \
    opihi_float *M2 = V2[0].vector[0].elementsInt; \
    for (i = 0; i < Nx; i++, out++, M1++, M2++) { \
  	*out = *M1 OP *M2; \
    } \
    break; \
  } \
  if ((M1->type != OPIHI_FLOAT) && (M2->type != OPIHI_FLOAT)) { \
    // output may be int, can use either as temp \
    opihi_float *M1 = V1[0].vector[0].elementsInt; \
    opihi_float *M2 = V2[0].vector[0].elementsInt; \
    for (i = 0; i < Nx; i++, out++, M1++, M2++) { \
  	*out = *M1 OP *M2; \
    } \
    break; \
  }

  switch (op[0]) {
  case '+': 
    VV_FUNC

    if ((M1->type == OPIHI_FLOAT) && (M2->type == OPIHI_FLOAT)) {
      // output must be float, can use either as a temp
      M1 = V1[0].vector[0].elements;
      M2 = V2[0].vector[0].elements;
      for (i = 0; i < Nx; i++, out++, M1++, M2++) {
	*out = *M1 + *M2;
      }
      break;
    }
    if ((M1->type != OPIHI_FLOAT) && (M2->type == OPIHI_FLOAT)) {
      // output must be float, can use M2 as a temp
      M1 = V1[0].vector[0].elementsInt;
      M2 = V2[0].vector[0].elements;
      for (i = 0; i < Nx; i++, out++, M1++, M2++) {
	*out = *M1 + *M2;
      }
      break;
    }
    if ((M1->type == OPIHI_FLOAT) && (M2->type != OPIHI_FLOAT)) {
      // output must be float, can use M1 as a temp
      M1 = V1[0].vector[0].elements;
      M2 = V2[0].vector[0].elementsInt;
      for (i = 0; i < Nx; i++, out++, M1++, M2++) {
	*out = *M1 + *M2;
      }
      break;
    }
    if ((M1->type != OPIHI_FLOAT) && (M2->type != OPIHI_FLOAT)) {
      // output may be int, can use either as temp
      M1 = V1[0].vector[0].elementsInt;
      M2 = V2[0].vector[0].elementsInt;
      for (i = 0; i < Nx; i++, out++, M1++, M2++) {
	*out = *M1 + *M2;
      }
      break;
    }
  case '-': 
    for (i = 0; i < Nx; i++, out++, M1++, M2++)
      *out = *M1 - *M2;
    break; 
  case '*': 
    for (i = 0; i < Nx; i++, out++, M1++, M2++)
      *out = *M1 * *M2;
    break; 
  case '/': 
    for (i = 0; i < Nx; i++, out++, M1++, M2++)
      *out = *M1 / *M2;
    break; 
  case '%': 
    for (i = 0; i < Nx; i++, out++, M1++, M2++)
      *out = (int) *M1 % (int) *M2;
    break; 
  case 0x5e: 
    for (i = 0; i < Nx; i++, out++, M1++, M2++) 
      *out = pow (*M1, *M2);
    break; 
  case '@': 
    for (i = 0; i < Nx; i++, out++, M1++, M2++) 
      *out = DEG_RAD*atan2 (*M1, *M2);
    break; 
  case 'D': 
    for (i = 0; i < Nx; i++, out++, M1++, M2++) 
      *out = MIN (*M1, *M2);
    break; 
  case 'U': 
    for (i = 0; i < Nx; i++, out++, M1++, M2++) 
      *out = MAX (*M1, *M2);
    break; 
  case '<': 
    for (i = 0; i < Nx; i++, out++, M1++, M2++) 
      *out = (*M1 < *M2) ? 1 : 0;
    break; 
  case '>': 
    for (i = 0; i < Nx; i++, out++, M1++, M2++) 
      *out = (*M1 > *M2) ? 1 : 0;
    break; 
  case '&': 
    for (i = 0; i < Nx; i++, out++, M1++, M2++) 
      *out = ((int)*M1 & (int)*M2);
    break; 
  case '|': 
    for (i = 0; i < Nx; i++, out++, M1++, M2++) 
      *out = ((int)*M1 | (int)*M2);
    break; 
  case 'E': 
    for (i = 0; i < Nx; i++, out++, M1++, M2++) 
      *out = (*M1 == *M2) ? 1 : 0;
    break; 
  case 'N': 
    for (i = 0; i < Nx; i++, out++, M1++, M2++) 
      *out = (*M1 != *M2) ? 1 : 0;
    break; 
  case 'L': 
    for (i = 0; i < Nx; i++, out++, M1++, M2++) 
      *out = (*M1 <= *M2) ? 1 : 0;
    break; 
  case 'G': 
    for (i = 0; i < Nx; i++, out++, M1++, M2++) 
      *out = (*M1 >= *M2) ? 1 : 0;
    break; 
  case 'A': 
    for (i = 0; i < Nx; i++, out++, M1++, M2++) 
      *out = (*M1 && *M2) ? 1 : 0;
    break; 
  case 'O': 
    for (i = 0; i < Nx; i++, out++, M1++, M2++) 
      *out = (*M1 || *M2) ? 1 : 0;
    break; 
  default:
    sprintf (line, "error: op %c not defined!", op[0]);
    push_error (line);
    return (FALSE);
  }

  /** free up any temporary buffers: **/

  if (V1[0].type == 'v') {
    free (V1[0].vector[0].elements);
    free (V1[0].vector);
  }
  if (V2[0].type == 'v') {
    free (V2[0].vector[0].elements);
    free (V2[0].vector);
  }
  /* at the end, V1 and V2 are deleted only if they were temporary */

  clear_stack (V1);
  clear_stack (V2);
  return (TRUE);

}

int SV_binary (StackVar *OUT, StackVar *V1, StackVar *V2, char *op) {

  int i, Nx;
  float *out, *M1, *M2;
  char line[512]; // this is only used to report an error 
  
  Nx = V2[0].vector[0].Nelements;

  if (V2[0].type == 'v') { /** use V2 as temp buffer, but header of V1 **/
    OUT[0].vector = V2[0].vector;
    V2[0].type = 'V'; /* prevent it from being freed below */
  } else {  /* no spare temp buffer */
    OUT[0].vector = InitVector ();
    CopyVector (OUT[0].vector, V2[0].vector);
  }
  OUT[0].type = 'v';   /*** <<--- says this is a temporary matrix ***/

  M1  = V1[0].ptr;
  M2  = V2[0].ptr;
  out = OUT[0].ptr = (float *)OUT[0].vector[0].elements;

  switch (op[0]) { 
  case '+': 
    for (i = 0; i < Nx; i++, out++, M2++)
      *out = *M1 + *M2;
    break; 
  case '-': 
    for (i = 0; i < Nx; i++, out++, M2++)
      *out = *M1 - *M2;
    break; 
  case '*': 
    for (i = 0; i < Nx; i++, out++, M2++)
      *out = *M1 * *M2;
    break; 
  case '/': 
    for (i = 0; i < Nx; i++, out++, M2++)
      *out = *M1 / *M2;
    break; 
  case '%': 
    for (i = 0; i < Nx; i++, out++, M2++)
      *out = (int) *M1 % (int) *M2;
    break; 
  case 0x5e: 
    for (i = 0; i < Nx; i++, out++, M2++) 
      *out = pow (*M1, *M2);
    break; 
  case '@': 
    for (i = 0; i < Nx; i++, out++, M2++) 
      *out = DEG_RAD*atan2 (*M1, *M2);
    break; 
  case 'D': 
    for (i = 0; i < Nx; i++, out++, M2++) 
      *out = MIN (*M1, *M2);
    break; 
  case 'U': 
    for (i = 0; i < Nx; i++, out++, M2++) 
      *out = MAX (*M1, *M2);
    break; 
  case '<': 
    for (i = 0; i < Nx; i++, out++, M2++) 
      *out = (*M1 < *M2) ? 1 : 0;
    break; 
  case '>': 
    for (i = 0; i < Nx; i++, out++, M2++) 
      *out = (*M1 > *M2) ? 1 : 0;
    break; 
  case '&': 
    for (i = 0; i < Nx; i++, out++, M2++) 
      *out = ((int)*M1 & (int)*M2);
    break; 
  case '|': 
    for (i = 0; i < Nx; i++, out++, M2++) 
      *out = ((int)*M1 | (int)*M2);
    break; 
  case 'E': 
    for (i = 0; i < Nx; i++, out++, M2++) 
      *out = (*M1 == *M2) ? 1 : 0;
    break; 
  case 'N': 
    for (i = 0; i < Nx; i++, out++, M2++) 
      *out = (*M1 != *M2) ? 1 : 0;
    break; 
  case 'L': 
    for (i = 0; i < Nx; i++, out++, M2++) 
      *out = (*M1 <= *M2) ? 1 : 0;
    break; 
  case 'G': 
    for (i = 0; i < Nx; i++, out++, M2++) 
      *out = (*M1 >= *M2) ? 1 : 0;
    break; 
  case 'A': 
    for (i = 0; i < Nx; i++, out++, M2++) 
      *out = (*M1 && *M2) ? 1 : 0;
    break; 
  case 'O': 
    for (i = 0; i < Nx; i++, out++, M2++) 
      *out = (*M1 || *M2) ? 1 : 0;
    break; 
  default:
    sprintf (line, "error: op %c not defined!", op[0]);
    push_error (line);
    return (FALSE);
  }

  /** free up any temporary buffers: **/
  if (V2[0].type == 'v') {
    free (V2[0].vector[0].elements);
    free (V2[0].vector);
  }

  clear_stack (V1);
  clear_stack (V2);

  /* at the end, V1 and V2 are deleted only if they were temporary */
  return (TRUE);

}

int VS_binary (StackVar *OUT, StackVar *V1, StackVar *V2, char *op) {

  int i, Nx;
  float *out, *M1, *M2;
  char line[512]; // this is only used to report an error 
  
  Nx = V1[0].vector[0].Nelements;

  if (V1[0].type == 'v') { /** use V1 as temp buffer **/
    OUT[0].vector = V1[0].vector;
    V1[0].type = 'V'; /* prevent it from being freed below */
  } else {  /* no spare temp buffer */
    OUT[0].vector = InitVector ();
    CopyVector (OUT[0].vector, V1[0].vector);
  }
  OUT[0].type = 'v';   /*** <<--- says this is a temporary matrix ***/

  M1  = V1[0].ptr;
  M2  = V2[0].ptr;
  out = OUT[0].ptr = (float *)OUT[0].vector[0].elements;

  switch (op[0]) { 
  case '+': 
    for (i = 0; i < Nx; i++, out++, M1++)
      *out = *M1 + *M2;
    break; 
  case '-': 
    for (i = 0; i < Nx; i++, out++, M1++)
      *out = *M1 - *M2;
    break; 
  case '*': 
    for (i = 0; i < Nx; i++, out++, M1++)
      *out = *M1 * *M2;
    break; 
  case '/': 
    for (i = 0; i < Nx; i++, out++, M1++)
      *out = *M1 / *M2;
    break; 
  case '%': 
    for (i = 0; i < Nx; i++, out++, M1++)
      *out = (int) *M1 % (int) *M2;
    break; 
  case 0x5e: 
    for (i = 0; i < Nx; i++, out++, M1++) 
      *out = pow (*M1, *M2);
    break; 
  case '@': 
    for (i = 0; i < Nx; i++, out++, M1++) 
      *out = DEG_RAD*atan2 (*M1, *M2);
    break; 
  case 'D': 
    for (i = 0; i < Nx; i++, out++, M1++) 
      *out = MIN (*M1, *M2);
    break; 
  case 'U': 
    for (i = 0; i < Nx; i++, out++, M1++) 
      *out = MAX (*M1, *M2);
    break; 
  case '<': 
    for (i = 0; i < Nx; i++, out++, M1++) 
      *out = (*M1 < *M2) ? 1 : 0;
    break; 
  case '>': 
    for (i = 0; i < Nx; i++, out++, M1++) 
      *out = (*M1 > *M2) ? 1 : 0;
    break; 
  case '&': 
    for (i = 0; i < Nx; i++, out++, M1++) 
      *out = ((int)*M1 & (int)*M2);
    break; 
  case '|': 
    for (i = 0; i < Nx; i++, out++, M1++) 
      *out = ((int)*M1 | (int)*M2);
    break; 
  case 'E': 
    for (i = 0; i < Nx; i++, out++, M1++) 
      *out = (*M1 == *M2) ? 1 : 0;
    break; 
  case 'N': 
    for (i = 0; i < Nx; i++, out++, M1++) 
      *out = (*M1 != *M2) ? 1 : 0;
    break; 
  case 'L': 
    for (i = 0; i < Nx; i++, out++, M1++) 
      *out = (*M1 <= *M2) ? 1 : 0;
    break; 
  case 'G': 
    for (i = 0; i < Nx; i++, out++, M1++) 
      *out = (*M1 >= *M2) ? 1 : 0;
    break; 
  case 'A': 
    for (i = 0; i < Nx; i++, out++, M1++) 
      *out = (*M1 && *M2) ? 1 : 0;
    break; 
  case 'O': 
    for (i = 0; i < Nx; i++, out++, M1++) 
      *out = (*M1 || *M2) ? 1 : 0;
    break; 
  default:
    sprintf (line, "error: op %c not defined!", op[0]);
    push_error (line);
    return (FALSE);
  }

  /** free up any temporary buffers: **/

  if (V1[0].type == 'v') {
    free (V1[0].vector[0].elements);
    free (V1[0].vector);
  }

  clear_stack (V1);
  clear_stack (V2);

  /* at the end, V1 and V2 are deleted only if they were temporary */
  return (TRUE);

}

int MV_binary (StackVar *OUT, StackVar *V1, StackVar *V2, char *op) {

  int i, j, Nx, Ny;
  float *out, *M1, *M2;
  char line[512]; // this is only used to report an error 
 
  Nx = V1[0].buffer[0].matrix.Naxis[0];
  Ny = V1[0].buffer[0].matrix.Naxis[1];
  if (Ny != V2[0].vector[0].Nelements) {
    push_error ("dimension mismatch");
    return (FALSE);
  }

  /* if possible, use V1 as temp buffer, otherwise create new one */
  if (V1[0].type == 'm') {  
    OUT[0].buffer = V1[0].buffer;
    V1[0].type = 'M'; /* prevent it from being freed below */
  } else {  
    /* do buffer.matrix.buffer and buffer.header.buffer get correctly zeroed? */
    OUT[0].buffer = InitBuffer ();
    CopyBuffer (OUT[0].buffer, V1[0].buffer);
  }
  OUT[0].type = 'm'; /*** <<--- says this is a temporary matrix ***/

  M1  = V1[0].ptr;
  M2  = V2[0].ptr;
  out = OUT[0].ptr = (float *)OUT[0].buffer[0].matrix.buffer;

  switch (op[0]) { 
  case '+': 
    for (i = 0; i < Ny; i++, M2++) {
      for (j = 0; j < Nx; j++, out++, M1++)
        *out = *M1 + *M2;
    }
    break; 
  case '-': 
    for (i = 0; i < Ny; i++, M2++) {
      for (j = 0; j < Nx; j++, out++, M1++)
        *out = *M1 - *M2;
    }
    break; 
  case '*': 
    for (i = 0; i < Ny; i++, M2++) {
      for (j = 0; j < Nx; j++, out++, M1++)
        *out = *M1 * *M2;
    }
    break; 
  case '/': 
    for (i = 0; i < Ny; i++, M2++) {
      for (j = 0; j < Nx; j++, out++, M1++)
        *out = *M1 / *M2;
    }
    break; 
  case '%': 
    for (i = 0; i < Ny; i++, M2++) {
      for (j = 0; j < Nx; j++, out++, M1++)
        *out = (int) *M1 % (int) *M2;
    }
    break; 
  case 0x5e: 
    for (i = 0; i < Ny; i++, M2++) {
      for (j = 0; j < Nx; j++, out++, M1++) 
        *out = pow (*M1, *M2);
    }
    break; 
  case '@': 
    for (i = 0; i < Ny; i++, M2++) {
      for (j = 0; j < Nx; j++, out++, M1++) 
        *out = DEG_RAD*atan2 (*M1, *M2);
    }
    break; 
  case 'D': 
    for (i = 0; i < Ny; i++, M2++) {
      for (j = 0; j < Nx; j++, out++, M1++) 
        *out = MIN (*M1, *M2);
    }
    break; 
  case 'U': 
    for (i = 0; i < Ny; i++, M2++) {
      for (j = 0; j < Nx; j++, out++, M1++) 
        *out = MAX (*M1, *M2);
    }
    break; 
  case '<': 
    for (i = 0; i < Ny; i++, M2++) {
      for (j = 0; j < Nx; j++, out++, M1++) 
        *out = (*M1 < *M2) ? 1 : 0;
    }
    break; 
  case '>': 
    for (i = 0; i < Ny; i++, M2++) {
      for (j = 0; j < Nx; j++, out++, M1++) 
        *out = (*M1 > *M2) ? 1 : 0;
    }
    break; 
  case '&': 
    for (i = 0; i < Ny; i++, M2++) {
      for (j = 0; j < Nx; j++, out++, M1++) 
        *out = ((int)*M1 & (int)*M2);
    }
    break; 
  case '|': 
    for (i = 0; i < Ny; i++, M2++) {
      for (j = 0; j < Nx; j++, out++, M1++) 
        *out = ((int)*M1 | (int)*M2);
    }
    break; 
  case 'E': 
    for (i = 0; i < Ny; i++, M2++) {
      for (j = 0; j < Nx; j++, out++, M1++) 
        *out = (*M1 == *M2) ? 1 : 0;
    }
    break; 
  case 'N': 
    for (i = 0; i < Ny; i++, M2++) {
      for (j = 0; j < Nx; j++, out++, M1++) 
        *out = (*M1 != *M2) ? 1 : 0;
    }
    break; 
  case 'L': 
    for (i = 0; i < Ny; i++, M2++) {
      for (j = 0; j < Nx; j++, out++, M1++) 
        *out = (*M1 <= *M2) ? 1 : 0;
    }
    break; 
  case 'G': 
    for (i = 0; i < Ny; i++, M2++) {
      for (j = 0; j < Nx; j++, out++, M1++) 
        *out = (*M1 >= *M2) ? 1 : 0;
    }
    break; 
  case 'A': 
    for (i = 0; i < Ny; i++, M2++) {
      for (j = 0; j < Nx; j++, out++, M1++) 
        *out = (*M1 && *M2) ? 1 : 0;
    }
    break; 
  case 'O': 
    for (i = 0; i < Ny; i++, M2++) {
      for (j = 0; j < Nx; j++, out++, M1++) 
        *out = (*M1 || *M2) ? 1 : 0;
    }
    break; 
  default:
    sprintf (line, "error: op %c not defined!", op[0]);
    push_error (line);
    return (FALSE);
  }

  /** free up any temporary buffers: **/

  if (V1[0].type == 'm') {
    free (V1[0].buffer[0].header.buffer);
    free (V1[0].buffer[0].matrix.buffer);
    free (V1[0].buffer);
  }
  if (V2[0].type == 'v') {
    free (V2[0].vector[0].elements);
    free (V2[0].vector);
  }

  clear_stack (V1);
  clear_stack (V2);

  /* at the end, V1 and V2 are deleted only if they were temporary */
  return (TRUE);

}


int VM_binary (StackVar *OUT, StackVar *V1, StackVar *V2, char *op) {

  int i, j, Nx, Ny;
  float *out, *M1, *M2;
  char line[512]; // this is only used to report an error 
  
  Nx = V2[0].buffer[0].matrix.Naxis[0];
  Ny = V2[0].buffer[0].matrix.Naxis[1];
  if (Nx != V1[0].vector[0].Nelements) {
    return (FALSE);
  }

  /* if possible, use V2 as temp buffer, otherwise create new one */
  if (V2[0].type == 'm') {
    OUT[0].buffer = V2[0].buffer;
    V2[0].type = 'M'; /* prevent it from being freed below */
  } else {  /* no spare temp buffer */
    OUT[0].buffer = InitBuffer ();
    CopyBuffer (OUT[0].buffer, V2[0].buffer);
  }
  OUT[0].type = 'm'; /*** <<--- says this is a temporary matrix ***/

  M1  = V1[0].ptr;
  M2  = V2[0].ptr;
  out = OUT[0].ptr = (float *)OUT[0].buffer[0].matrix.buffer;

  switch (op[0]) { 
  case '+': 
    for (i = 0; i < Ny; i++) {
      M1 = V1[0].ptr;
      for (j = 0; j < Nx; j++, out++, M1++, M2++)
        *out = *M1 + *M2;
    }
    break; 
  case '-': 
    for (i = 0; i < Ny; i++) {
      M1 = V1[0].ptr;
      for (j = 0; j < Nx; j++, out++, M1++, M2++)
        *out = *M1 - *M2;
    }
    break; 
  case '*': 
    for (i = 0; i < Ny; i++) {
      M1 = V1[0].ptr;
      for (j = 0; j < Nx; j++, out++, M1++, M2++)
        *out = *M1 * *M2;
    }
    break; 
  case '/': 
    for (i = 0; i < Ny; i++) {
      M1 = V1[0].ptr;
      for (j = 0; j < Nx; j++, out++, M1++, M2++)
        *out = *M1 / *M2;
    }
    break; 
  case '%': 
    for (i = 0; i < Ny; i++) {
      M1 = V1[0].ptr;
      for (j = 0; j < Nx; j++, out++, M1++, M2++)
        *out = (int) *M1 % (int) *M2;
    }
    break; 
  case 0x5e: 
    for (i = 0; i < Ny; i++) {
      M1 = V1[0].ptr;
      for (j = 0; j < Nx; j++, out++, M1++, M2++) 
        *out = pow (*M1, *M2);
    }
    break; 
  case '@': 
    for (i = 0; i < Ny; i++) {
      M1 = V1[0].ptr;
      for (j = 0; j < Nx; j++, out++, M1++, M2++) 
        *out = DEG_RAD*atan2 (*M1, *M2);
    }
    break; 
  case 'D': 
    for (i = 0; i < Ny; i++) {
      M1 = V1[0].ptr;
      for (j = 0; j < Nx; j++, out++, M1++, M2++) 
        *out = MIN (*M1, *M2);
    }
    break; 
  case 'U': 
    for (i = 0; i < Ny; i++) {
      M1 = V1[0].ptr;
      for (j = 0; j < Nx; j++, out++, M1++, M2++) 
        *out = MAX (*M1, *M2);
    }
    break; 
  case '<': 
    for (i = 0; i < Ny; i++) {
      M1 = V1[0].ptr;
      for (j = 0; j < Nx; j++, out++, M1++, M2++) 
        *out = (*M1 < *M2) ? 1 : 0;
    }
    break; 
  case '>': 
    for (i = 0; i < Ny; i++) {
      M1 = V1[0].ptr;
      for (j = 0; j < Nx; j++, out++, M1++, M2++) 
        *out = (*M1 > *M2) ? 1 : 0;
    }
    break; 
  case '&': 
    for (i = 0; i < Ny; i++) {
      M1 = V1[0].ptr;
      for (j = 0; j < Nx; j++, out++, M1++, M2++) 
        *out = ((int)*M1 & (int)*M2);
    }
    break; 
  case '|': 
    for (i = 0; i < Ny; i++) {
      M1 = V1[0].ptr;
      for (j = 0; j < Nx; j++, out++, M1++, M2++) 
        *out = ((int)*M1 | (int)*M2);
    }
    break; 
  case 'E': 
    for (i = 0; i < Ny; i++) {
      M1 = V1[0].ptr;
      for (j = 0; j < Nx; j++, out++, M1++, M2++) 
        *out = (*M1 == *M2) ? 1 : 0;
    }
    break; 
  case 'N': 
    for (i = 0; i < Ny; i++) {
      M1 = V1[0].ptr;
      for (j = 0; j < Nx; j++, out++, M1++, M2++) 
        *out = (*M1 != *M2) ? 1 : 0;
    }
    break; 
  case 'L': 
    for (i = 0; i < Ny; i++) {
      M1 = V1[0].ptr;
      for (j = 0; j < Nx; j++, out++, M1++, M2++) 
        *out = (*M1 <= *M2) ? 1 : 0;
    }
    break; 
  case 'G': 
    for (i = 0; i < Ny; i++) {
      M1 = V1[0].ptr;
      for (j = 0; j < Nx; j++, out++, M1++, M2++) 
        *out = (*M1 >= *M2) ? 1 : 0;
    }
    break; 
  case 'A': 
    for (i = 0; i < Ny; i++) {
      M1 = V1[0].ptr;
      for (j = 0; j < Nx; j++, out++, M1++, M2++) 
        *out = (*M1 && *M2) ? 1 : 0;
    }
    break; 
  case 'O': 
    for (i = 0; i < Ny; i++) {
      M1 = V1[0].ptr;
      for (j = 0; j < Nx; j++, out++, M1++, M2++) 
        *out = (*M1 || *M2) ? 1 : 0;
    }
    break; 
  default:
    sprintf (line, "error: op %c not defined!", op[0]);
    push_error (line);
    return (FALSE);
  }

  /** free up any temporary buffers: **/

  if (V1[0].type == 'v') {
    free (V1[0].vector[0].elements);
    free (V1[0].vector);
  }
  if (V2[0].type == 'm') {
    free (V2[0].buffer[0].header.buffer);
    free (V2[0].buffer[0].matrix.buffer);
    free (V2[0].buffer);
  }

  clear_stack (V1);
  clear_stack (V2);

  /* at the end, V1 and V2 are deleted only if they were temporary */
  return (TRUE);

}

int MM_binary (StackVar *OUT, StackVar *V1, StackVar *V2, char *op) {

  int i, Nx, Ny;
  float *out, *M1, *M2;
  char line[512]; // this is only used to report an error 
  
  Nx = V1[0].buffer[0].matrix.Naxis[0];
  Ny = V1[0].buffer[0].matrix.Naxis[1];

  if (V1[0].type == 'm') {  /** use V1 as temp buffer **/
    OUT[0].buffer = V1[0].buffer;
    V1[0].type = 'M'; /* prevent it from being freed below */
  } else {
    if (V2[0].type == 'm') { /** use V2 as temp buffer, but header of V1 **/
      OUT[0].buffer = V2[0].buffer;
      V2[0].type = 'M'; /* prevent it from being freed below */
    } else {  /* no spare temp buffer */
      OUT[0].buffer = InitBuffer ();
      CopyBuffer (OUT[0].buffer, V1[0].buffer);
    }
  }
  OUT[0].type = 'm'; /*** <<--- says this is a temporary matrix ***/

  M1  = V1[0].ptr;
  M2  = V2[0].ptr;
  out = OUT[0].ptr = (float *)OUT[0].buffer[0].matrix.buffer;

  switch (op[0]) { 
  case '+': 
    for (i = 0; i < Nx*Ny; i++, out++, M1++, M2++)
      *out = *M1 + *M2;
    break; 
  case '-': 
    for (i = 0; i < Nx*Ny; i++, out++, M1++, M2++)
      *out = *M1 - *M2;
    break; 
  case '*': 
    for (i = 0; i < Nx*Ny; i++, out++, M1++, M2++)
      *out = *M1 * *M2;
    break; 
  case '/': 
    for (i = 0; i < Nx*Ny; i++, out++, M1++, M2++)
      *out = *M1 / *M2;
    break; 
  case '%': 
    for (i = 0; i < Nx*Ny; i++, out++, M1++, M2++)
      *out = (int) *M1 % (int) *M2;
    break; 
  case 0x5e: 
    for (i = 0; i < Nx*Ny; i++, out++, M1++, M2++) 
      *out = pow (*M1, *M2);
    break; 
  case '@': 
    for (i = 0; i < Nx*Ny; i++, out++, M1++, M2++) 
      *out = DEG_RAD*atan2 (*M1, *M2);
    break; 
  case 'D': 
    for (i = 0; i < Nx*Ny; i++, out++, M1++, M2++) 
      *out = MIN (*M1, *M2);
    break; 
  case 'U': 
    for (i = 0; i < Nx*Ny; i++, out++, M1++, M2++) 
      *out = MAX (*M1, *M2);
    break; 
  case '<': 
    for (i = 0; i < Nx*Ny; i++, out++, M1++, M2++) 
      *out = (*M1 < *M2) ? 1 : 0;
    break; 
  case '>': 
    for (i = 0; i < Nx*Ny; i++, out++, M1++, M2++) 
      *out = (*M1 > *M2) ? 1 : 0;
    break; 
  case '&': 
    for (i = 0; i < Nx*Ny; i++, out++, M1++, M2++) 
      *out = ((int)*M1 & (int)*M2);
    break; 
  case '|': 
    for (i = 0; i < Nx*Ny; i++, out++, M1++, M2++) 
      *out = ((int)*M1 | (int)*M2);
    break; 
  case 'E': 
    for (i = 0; i < Nx*Ny; i++, out++, M1++, M2++) 
      *out = (*M1 == *M2) ? 1 : 0;
    break; 
  case 'N': 
    for (i = 0; i < Nx*Ny; i++, out++, M1++, M2++) 
      *out = (*M1 != *M2) ? 1 : 0;
    break; 
  case 'L': 
    for (i = 0; i < Nx*Ny; i++, out++, M1++, M2++) 
      *out = (*M1 <= *M2) ? 1 : 0;
    break; 
  case 'G': 
    for (i = 0; i < Nx*Ny; i++, out++, M1++, M2++) 
      *out = (*M1 >= *M2) ? 1 : 0;
    break; 
  case 'A': 
    for (i = 0; i < Nx*Ny; i++, out++, M1++, M2++) 
      *out = (*M1 && *M2) ? 1 : 0;
    break; 
  case 'O': 
    for (i = 0; i < Nx*Ny; i++, out++, M1++, M2++) 
      *out = (*M1 || *M2) ? 1 : 0;
    break; 
  default:
    sprintf (line, "error: op %c not defined!", op[0]);
    push_error (line);
    return (FALSE);
  }

  /** free up any temporary buffers: **/

  if (V1[0].type == 'm') {
    free (V1[0].buffer[0].header.buffer);
    free (V1[0].buffer[0].matrix.buffer);
    free (V1[0].buffer);
  }
  if (V2[0].type == 'm') {
    free (V2[0].buffer[0].header.buffer);
    free (V2[0].buffer[0].matrix.buffer);
    free (V2[0].buffer);
  }

  clear_stack (V1);
  clear_stack (V2);

  /* at the end, V1 and V2 are deleted only if they were temporary */
  return (TRUE);

}


int MS_binary (StackVar *OUT, StackVar *V1, StackVar *V2, char *op) {

  int i, Nx, Ny;
  float *out, *M1, *M2;
  char line[512]; // this is only used to report an error 
  
  Nx = V1[0].buffer[0].matrix.Naxis[0];
  Ny = V1[0].buffer[0].matrix.Naxis[1];

  /* if possible, use V1 as temp buffer, otherwise create new one */
  if (V1[0].type == 'm') {
    OUT[0].buffer = V1[0].buffer;
    V1[0].type = 'M'; /* prevent it from being freed below */
  } else {
    OUT[0].buffer = InitBuffer ();
    CopyBuffer (OUT[0].buffer, V1[0].buffer);
  }
  OUT[0].type = 'm';      /*** <<--- says this is a temporary matrix ***/

  M1  = V1[0].ptr;
  M2  = V2[0].ptr;
  out = OUT[0].ptr = (float *)OUT[0].buffer[0].matrix.buffer;

  switch (op[0]) { 
  case '+': 
    for (i = 0; i < Nx*Ny; i++, out++, M1++)
      *out = *M1 + *M2;
    break; 
  case '-': 
    for (i = 0; i < Nx*Ny; i++, out++, M1++)
      *out = *M1 - *M2;
    break; 
  case '*': 
    for (i = 0; i < Nx*Ny; i++, out++, M1++)
      *out = *M1 * *M2;
    break; 
  case '/': 
    for (i = 0; i < Nx*Ny; i++, out++, M1++)
      *out = *M1 / *M2;
    break; 
  case '%': 
    for (i = 0; i < Nx*Ny; i++, out++, M1++)
      *out = (int) *M1 % (int) *M2;
    break; 
  case 0x5e: 
    for (i = 0; i < Nx*Ny; i++, out++, M1++) 
      *out = pow (*M1, *M2);
    break; 
  case '@': 
    for (i = 0; i < Nx*Ny; i++, out++, M1++) 
      *out = DEG_RAD*atan2 (*M1, *M2);
    break; 
  case 'D': 
    for (i = 0; i < Nx*Ny; i++, out++, M1++) 
      *out = MIN (*M1, *M2);
    break; 
  case 'U': 
    for (i = 0; i < Nx*Ny; i++, out++, M1++) 
      *out = MAX (*M1, *M2);
    break; 
  case '<': 
    for (i = 0; i < Nx*Ny; i++, out++, M1++) 
      *out = (*M1 < *M2) ? 1 : 0;
    break; 
  case '>': 
    for (i = 0; i < Nx*Ny; i++, out++, M1++) 
      *out = (*M1 > *M2) ? 1 : 0;
    break; 
  case '&': 
    for (i = 0; i < Nx*Ny; i++, out++, M1++) 
      *out = ((int)*M1 & (int)*M2);
    break; 
  case '|': 
    for (i = 0; i < Nx*Ny; i++, out++, M1++) 
      *out = ((int)*M1 | (int)*M2);
    break; 
  case 'E': 
    for (i = 0; i < Nx*Ny; i++, out++, M1++) 
      *out = (*M1 == *M2) ? 1 : 0;
    break; 
  case 'N': 
    for (i = 0; i < Nx*Ny; i++, out++, M1++) 
      *out = (*M1 != *M2) ? 1 : 0;
    break; 
  case 'L': 
    for (i = 0; i < Nx*Ny; i++, out++, M1++) 
      *out = (*M1 <= *M2) ? 1 : 0;
    break; 
  case 'G': 
    for (i = 0; i < Nx*Ny; i++, out++, M1++) 
      *out = (*M1 >= *M2) ? 1 : 0;
    break; 
  case 'A': 
    for (i = 0; i < Nx*Ny; i++, out++, M1++) 
      *out = (*M1 && *M2) ? 1 : 0;
    break; 
  case 'O': 
    for (i = 0; i < Nx*Ny; i++, out++, M1++) 
      *out = (*M1 || *M2) ? 1 : 0;
    break; 
  default:
    sprintf (line, "error: op %c not defined!", op[0]);
    push_error (line);
    return (FALSE);
  }

  if (V1[0].type == 'm') {
    free (V1[0].buffer[0].header.buffer);
    free (V1[0].buffer[0].matrix.buffer);
    free (V1[0].buffer);
  }
  clear_stack (V1);
  clear_stack (V2);

  return (TRUE);

}


int SM_binary (StackVar *OUT, StackVar *V1, StackVar *V2, char *op) {

  int i, Nx, Ny;
  float *out, *M1, *M2;
  char line[512]; // this is only used to report an error 
  
  Nx = V2[0].buffer[0].matrix.Naxis[0];
  Ny = V2[0].buffer[0].matrix.Naxis[1];

  if (V2[0].type == 'm') {  /* V2[0] is NOT temporary, we can't use it for storage */
    OUT[0].buffer = V2[0].buffer;
    V2[0].type = 'M'; /* prevent it from being freed below */
  } else {
    OUT[0].buffer = InitBuffer ();
    CopyBuffer (OUT[0].buffer, V2[0].buffer);
  }
  OUT[0].type = 'm'; /*** <<--- says this is a temporary matrix ***/

  M1  = V1[0].ptr;
  M2  = V2[0].ptr;
  out = OUT[0].ptr = (float *)OUT[0].buffer[0].matrix.buffer;

  switch (op[0]) { 
  case '+': 
    for (i = 0; i < Nx*Ny; i++, out++, M2++)
      *out = *M1 + *M2;
    break; 
  case '-': 
    for (i = 0; i < Nx*Ny; i++, out++, M2++)
      *out = *M1 - *M2;
    break; 
  case '*': 
    for (i = 0; i < Nx*Ny; i++, out++, M2++)
      *out = *M1 * *M2;
    break; 
  case '/': 
    for (i = 0; i < Nx*Ny; i++, out++, M2++)
      *out = *M1 / *M2;
    break; 
  case '%': 
    for (i = 0; i < Nx*Ny; i++, out++, M2++)
      *out = (int) *M1 % (int) *M2;
    break; 
  case 0x5e: 
    for (i = 0; i < Nx*Ny; i++, out++, M2++)
      *out = pow (*M1, *M2);
    break; 
  case '@': 
    for (i = 0; i < Nx*Ny; i++, out++, M2++) 
      *out = DEG_RAD*atan2 (*M1, *M2);
    break; 
  case 'D': 
    for (i = 0; i < Nx*Ny; i++, out++, M2++) 
      *out = MIN (*M1, *M2);
    break; 
  case 'U': 
    for (i = 0; i < Nx*Ny; i++, out++, M2++) 
      *out = MAX (*M1, *M2);
    break; 
  case '<': 
    for (i = 0; i < Nx*Ny; i++, out++, M2++) 
      *out = (*M1 < *M2) ? 1 : 0;
    break; 
  case '>': 
    for (i = 0; i < Nx*Ny; i++, out++, M2++) 
      *out = (*M1 > *M2) ? 1 : 0;
    break; 
  case '&': 
    for (i = 0; i < Nx*Ny; i++, out++, M2++) 
      *out = ((int)*M1 & (int)*M2);
    break; 
  case '|': 
    for (i = 0; i < Nx*Ny; i++, out++, M2++) 
      *out = ((int)*M1 | (int)*M2);
    break; 
  case 'E': 
    for (i = 0; i < Nx*Ny; i++, out++, M2++) 
      *out = (*M1 == *M2) ? 1 : 0;
    break; 
  case 'N': 
    for (i = 0; i < Nx*Ny; i++, out++, M2++) 
      *out = (*M1 != *M2) ? 1 : 0;
    break; 
  case 'L': 
    for (i = 0; i < Nx*Ny; i++, out++, M2++) 
      *out = (*M1 <= *M2) ? 1 : 0;
    break; 
  case 'G': 
    for (i = 0; i < Nx*Ny; i++, out++, M2++) 
      *out = (*M1 >= *M2) ? 1 : 0;
    break; 
  case 'A': 
    for (i = 0; i < Nx*Ny; i++, out++, M2++) 
      *out = (*M1 && *M2) ? 1 : 0;
    break; 
  case 'O': 
    for (i = 0; i < Nx*Ny; i++, out++, M2++) 
      *out = (*M1 || *M2) ? 1 : 0;
    break; 
  default:
    sprintf (line, "error: op %c not defined!", op[0]);
    push_error (line);
    return (FALSE);
  }

  if (V2[0].type == 'm') {
    free (V2[0].buffer[0].header.buffer);
    free (V2[0].buffer[0].matrix.buffer);
    free (V2[0].buffer);
  }
  clear_stack (V1);
  clear_stack (V2);

  return (TRUE);

}


int SS_binary (StackVar *OUT, StackVar *V1, StackVar *V2, char *op) {

  float *M1, *M2, *out;
  char line[512]; // this is only used to report an error 

  M1  = V1[0].ptr;
  M2  = V2[0].ptr;
  OUT[0].ptr = V1[0].ptr;
  out = OUT[0].ptr;


  switch (op[0]) { 
  case '+': 
    *out = *M1 + *M2;
    break;    
  case '-': 
    *out = *M1 - *M2;
    break;    
  case '*': 
    *out = *M1 * *M2;
    break;    
  case '/': 
    *out = *M1 / *M2;
    break; 
  case '%': 
    *out = (int) *M1 % (int) *M2;
    break; 
  case 0x5e: 
    *out = pow (*M1, *M2);
    break; 
  case '@': 
    *out = DEG_RAD*atan2 (*M1, *M2);
    break; 
  case 'D': 
    *out = MIN (*M1, *M2);
    break; 
  case 'U': 
    *out = MAX (*M1, *M2);
    break; 
  case '<': 
    *out = (*M1 < *M2) ? 1 : 0;
    break; 
  case '>': 
    *out = (*M1 > *M2) ? 1 : 0;
    break; 
  case '&': 
    *out = ((int)*M1 & (int)*M2);
    break; 
  case '|': 
    *out = ((int)*M1 | (int)*M2);
    break; 
  case 'E': 
    *out = (*M1 == *M2) ? 1 : 0;
    break; 
  case 'N': 
    *out = (*M1 != *M2) ? 1 : 0;
    break; 
  case 'L': 
    *out = (*M1 <= *M2) ? 1 : 0;
    break; 
  case 'G': 
    *out = (*M1 >= *M2) ? 1 : 0;
    break; 
  case 'A': 
    *out = (*M1 && *M2) ? 1 : 0;
    break; 
  case 'O': 
    *out = (*M1 || *M2) ? 1 : 0;
    break; 
  default:
    sprintf (line, "error: op %c not defined!", op[0]);
    push_error (line);
    return (FALSE);
  }
  OUT[0].Float = *(OUT[0].ptr);
  OUT[0].type = 'S';

  clear_stack (V1);
  clear_stack (V2);

  return (TRUE);

}

int WW_binary (StackVar *OUT, StackVar *V1, StackVar *V2, char *op) {

  int value;
  char line[512]; // this is only used to report an error 

  /* only 'N' and 'E' are allowed for WW_binary operations. anything else is either a
   syntax error or is a string which looks like a math expression. */

  if ((op[0] != 'N') && (op[0] != 'E')) {
    sprintf (line, "error: op %c not defined for string operations!", op[0]);
    push_error (line);
    return (FALSE);
  }

  /* evaluate stack will only call WW_binary with one numerical value,
     and only in the case that the string did not parse to a number 
     thus: string == number -> false */

  if (!strncasecmp (&V1[0].type, "S", 1)) {
    value = (op[0] == 'N');
    goto escape;
  }
  if (!strncasecmp (&V2[0].type, "S", 1)) {
    value = (op[0] == 'N');
    goto escape;
  }

  switch (op[0]) { 
  case 'E': 
    value = strcmp (V1[0].name, V2[0].name) ? 0 : 1;
    break; 
  case 'N': 
    value = strcmp (V1[0].name, V2[0].name) ? 1 : 0;
    break; 
  default:
    sprintf (line, "error: op %c not defined for string operations!", op[0]);
    push_error (line);
    return (FALSE);
  }

escape:
  OUT[0].Float = value;
  OUT[0].type = 'S';
  OUT[0].ptr = &OUT[0].Float;

  clear_stack (V1);
  clear_stack (V2);
  return (TRUE);

}


int S_unary (StackVar *OUT, StackVar *V1, char *op) {

  float *out, *M1;
  
  out = OUT[0].ptr = V1[0].ptr;
  M1  = V1[0].ptr;

  if (!strcmp (op, "="))     {    }
  if (!strcmp (op, "abs"))   {    *out = fabs(*M1);           }
  if (!strcmp (op, "int"))   {    *out = (float)(int)(*M1);   }
  if (!strcmp (op, "exp"))   {    *out = exp (*M1);           }
  if (!strcmp (op, "ten"))   {    *out = pow (10.0,*M1);      }
  if (!strcmp (op, "log"))   {    *out = log10 (*M1);         }
  if (!strcmp (op, "ln"))    {    *out = log (*M1);           }
  if (!strcmp (op, "sqrt"))  {    *out = sqrt (*M1);          }
  if (!strcmp (op, "erf"))   {    *out = erf (*M1);           }

  if (!strcmp (op, "sinh"))  {    *out = sinh (*M1);          }
  if (!strcmp (op, "cosh"))  {    *out = cosh (*M1);          }
  if (!strcmp (op, "asinh")) {    *out = asinh (*M1);         }
  if (!strcmp (op, "acosh")) {    *out = acosh (*M1);         }
  if (!strcmp (op, "lgamma")) {   *out = lgamma (*M1);        }

  if (!strcmp (op, "sin"))   {    *out = sin (*M1);           }
  if (!strcmp (op, "cos"))   {    *out = cos (*M1);           }
  if (!strcmp (op, "tan"))   {    *out = tan (*M1);           }
  if (!strcmp (op, "dsin"))  {    *out = sin (*M1*RAD_DEG);   }
  if (!strcmp (op, "dcos"))  {    *out = cos (*M1*RAD_DEG);   }
  if (!strcmp (op, "dtan"))  {    *out = tan (*M1*RAD_DEG);   }
  if (!strcmp (op, "asin"))  {    *out = asin (*M1);          }
  if (!strcmp (op, "acos"))  {    *out = acos (*M1);          }
  if (!strcmp (op, "atan"))  {    *out = atan (*M1);          }
  if (!strcmp (op, "dasin")) {    *out = asin (*M1)*DEG_RAD;  }
  if (!strcmp (op, "dacos")) {    *out = acos (*M1)*DEG_RAD;  }
  if (!strcmp (op, "datan")) {    *out = atan (*M1)*DEG_RAD;  }
  if (!strcmp (op, "rnd"))   {    *out = drand48();           }
  if (!strcmp (op, "not"))   {    *out = !(*M1);              }
  if (!strcmp (op, "--"))    {    *out = - (*M1);             }
  if (!strcmp (op, "isinf")) {    *out = !finite(*M1);        }
  if (!strcmp (op, "isnan")) {    *out = isnan(*M1);          } 

  OUT[0].Float = *out;
  OUT[0].type = 'S';

  clear_stack (V1);
  return (TRUE);

}

int V_unary (StackVar *OUT, StackVar *V1, char *op) {

  int i, Nx;
  float *out, *M1;
  
  Nx = V1[0].vector[0].Nelements;
 
  if (V1[0].type == 'v') {  /** use V1 as temp buffer **/
    OUT[0].vector = V1[0].vector;
    V1[0].type = 'V'; /* prevent it from being freed below */
  } else {  /* no spare temp buffer */
    OUT[0].vector = InitVector ();
    CopyVector (OUT[0].vector, V1[0].vector);
  }
  OUT[0].type = 'v'; /*** <<--- says this is a temporary matrix ***/
  M1  = V1[0].ptr;
  out = OUT[0].ptr = (float *)OUT[0].vector[0].elements;

  if (!strcmp (op, "="))     { } /* already set equal */
  if (!strcmp (op, "abs"))   { for (i = 0; i < Nx; i++, out++, M1++) { *out = fabs(*M1);         }}
  if (!strcmp (op, "int"))   { for (i = 0; i < Nx; i++, out++, M1++) { *out = (float)(int)(*M1); }}
  if (!strcmp (op, "exp"))   { for (i = 0; i < Nx; i++, out++, M1++) { *out = exp(*M1);          }}
  if (!strcmp (op, "ten"))   { for (i = 0; i < Nx; i++, out++, M1++) { *out = pow(10.0,*M1);     }}
  if (!strcmp (op, "log"))   { for (i = 0; i < Nx; i++, out++, M1++) { *out = log10(*M1);        }}
  if (!strcmp (op, "ln"))    { for (i = 0; i < Nx; i++, out++, M1++) { *out = log(*M1);          }}
  if (!strcmp (op, "sqrt"))  { for (i = 0; i < Nx; i++, out++, M1++) { *out = sqrt(*M1);         }}
  if (!strcmp (op, "erf"))   { for (i = 0; i < Nx; i++, out++, M1++) { *out = erf(*M1);          }}

  if (!strcmp (op, "sinh"))   { for (i = 0; i < Nx; i++, out++, M1++) { *out = sinh(*M1);        }}
  if (!strcmp (op, "cosh"))   { for (i = 0; i < Nx; i++, out++, M1++) { *out = cosh(*M1);        }}
  if (!strcmp (op, "asinh"))  { for (i = 0; i < Nx; i++, out++, M1++) { *out = asinh(*M1);       }}
  if (!strcmp (op, "acosh"))  { for (i = 0; i < Nx; i++, out++, M1++) { *out = acosh(*M1);       }}
  if (!strcmp (op, "lgamma")) { for (i = 0; i < Nx; i++, out++, M1++) { *out = lgamma(*M1);      }}

  if (!strcmp (op, "sin"))   { for (i = 0; i < Nx; i++, out++, M1++) { *out = sin(*M1);          }}
  if (!strcmp (op, "cos"))   { for (i = 0; i < Nx; i++, out++, M1++) { *out = cos(*M1);          }}
  if (!strcmp (op, "tan"))   { for (i = 0; i < Nx; i++, out++, M1++) { *out = tan(*M1);          }}
  if (!strcmp (op, "dsin"))  { for (i = 0; i < Nx; i++, out++, M1++) { *out = sin(*M1*RAD_DEG);  }}
  if (!strcmp (op, "dcos"))  { for (i = 0; i < Nx; i++, out++, M1++) { *out = cos(*M1*RAD_DEG);  }}
  if (!strcmp (op, "dtan"))  { for (i = 0; i < Nx; i++, out++, M1++) { *out = tan(*M1*RAD_DEG);  }}
  if (!strcmp (op, "asin"))  { for (i = 0; i < Nx; i++, out++, M1++) { *out = asin(*M1);         }}
  if (!strcmp (op, "acos"))  { for (i = 0; i < Nx; i++, out++, M1++) { *out = acos(*M1);         }}
  if (!strcmp (op, "atan"))  { for (i = 0; i < Nx; i++, out++, M1++) { *out = atan(*M1);         }}
  if (!strcmp (op, "dasin")) { for (i = 0; i < Nx; i++, out++, M1++) { *out = asin(*M1)*DEG_RAD; }}
  if (!strcmp (op, "dacos")) { for (i = 0; i < Nx; i++, out++, M1++) { *out = acos(*M1)*DEG_RAD; }}
  if (!strcmp (op, "datan")) { for (i = 0; i < Nx; i++, out++, M1++) { *out = atan(*M1)*DEG_RAD; }}
  if (!strcmp (op, "rnd"))   { for (i = 0; i < Nx; i++, out++, M1++) { *out = drand48();         }}
  if (!strcmp (op, "ramp"))  { for (i = 0; i < Nx; i++, out++, M1++) { *out = i;                 }}
  if (!strcmp (op, "zero"))  { for (i = 0; i < Nx; i++, out++, M1++) { *out = 0;                 }}
  if (!strcmp (op, "not"))   { for (i = 0; i < Nx; i++, out++, M1++) { *out = !(*M1);            }}
  if (!strcmp (op, "--"))    { for (i = 0; i < Nx; i++, out++, M1++) { *out = -(*M1);            }}
  if (!strcmp (op, "isinf")) { for (i = 0; i < Nx; i++, out++, M1++) { *out = !finite(*M1);      }}
  if (!strcmp (op, "isnan")) { for (i = 0; i < Nx; i++, out++, M1++) { *out = isnan(*M1);        }}

  /* xramp and yramp only make sense in for matrices. for vectors, xramp = ramp, yramp = zero */
  if (!strcmp (op, "xramp")) { for (i = 0; i < Nx; i++, out++, M1++) { *out = i;                 }}
  if (!strcmp (op, "yramp")) { for (i = 0; i < Nx; i++, out++, M1++) { *out = 0;                 }}

  if (V1[0].type == 'v') {
    free (V1[0].vector[0].elements);
    free (V1[0].vector);
    V1[0].vector = NULL;
  }  

  clear_stack (V1);
  return (TRUE);

}

int M_unary (StackVar *OUT, StackVar *V1, char *op) {

  int i, j, Nx, Ny;
  float *out, *M1;
  
  Nx = V1[0].buffer[0].matrix.Naxis[0];
  Ny = V1[0].buffer[0].matrix.Naxis[1];

  if (V1[0].type == 'm') {
    OUT[0].buffer = V1[0].buffer;
    V1[0].type = 'M'; /* prevent it from being freed below */
  } else {
    OUT[0].buffer = InitBuffer ();
    CopyBuffer (OUT[0].buffer, V1[0].buffer);
  }
  OUT[0].type = 'm';      /*** <<--- says this is a temporary matrix ***/
  M1  = V1[0].ptr;
  out = OUT[0].ptr = (float *)OUT[0].buffer[0].matrix.buffer;

  if (!strcmp (op, "="))     { }
  if (!strcmp (op, "abs"))   { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = fabs(*M1);         }}
  if (!strcmp (op, "int"))   { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = (float)(int)(*M1); }}
  if (!strcmp (op, "exp"))   { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = exp(*M1);          }}
  if (!strcmp (op, "ten"))   { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = pow(10.0,*M1);     }}
  if (!strcmp (op, "log"))   { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = log10(*M1);        }}
  if (!strcmp (op, "ln"))    { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = log(*M1);          }}
  if (!strcmp (op, "sqrt"))  { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = sqrt(*M1);         }}
  if (!strcmp (op, "erf"))   { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = erf(*M1);          }}

  if (!strcmp (op, "sinh"))  { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = sinh(*M1);         }}
  if (!strcmp (op, "cosh"))  { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = cosh(*M1);         }}
  if (!strcmp (op, "asinh")) { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = asinh(*M1);        }}
  if (!strcmp (op, "acosh")) { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = acosh(*M1);        }}
  if (!strcmp (op, "lgamma")) { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = lgamma(*M1);      }}

  if (!strcmp (op, "sin"))   { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = sin(*M1);          }}
  if (!strcmp (op, "cos"))   { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = cos(*M1);          }}
  if (!strcmp (op, "tan"))   { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = tan(*M1);          }}
  if (!strcmp (op, "dsin"))  { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = sin(*M1*RAD_DEG);  }}
  if (!strcmp (op, "dcos"))  { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = cos(*M1*RAD_DEG);  }}
  if (!strcmp (op, "dtan"))  { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = tan(*M1*RAD_DEG);  }}
  if (!strcmp (op, "asin"))  { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = asin(*M1);         }}
  if (!strcmp (op, "acos"))  { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = acos(*M1);         }}
  if (!strcmp (op, "atan"))  { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = atan(*M1);         }}
  if (!strcmp (op, "dasin")) { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = asin(*M1)*DEG_RAD; }}
  if (!strcmp (op, "dacos")) { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = acos(*M1)*DEG_RAD; }}
  if (!strcmp (op, "datan")) { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = atan(*M1)*DEG_RAD; }}
  if (!strcmp (op, "not"))   { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = !(*M1);            }}
  if (!strcmp (op, "--"))    { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = -(*M1);            }}
  if (!strcmp (op, "rnd"))   { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = drand48();         }}
  if (!strcmp (op, "ramp"))  { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = i;                 }}
  if (!strcmp (op, "zero"))  { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = 0;                 }}
  if (!strcmp (op, "isinf")) { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = !finite(*M1);      }}
  if (!strcmp (op, "isnan")) { for (i = 0; i < Nx*Ny; i++, out++, M1++) { *out = isnan(*M1);        }}

  /* xrm and yrm only make sense in for matrices. see special meaning for vectors */
  if (!strcmp (op, "xramp")) {
    for (j = 0; j < Ny; j++) {
      for (i = 0; i < Nx; i++, out++, M1++) {
        *out = i;
      }
    }
  }
  if (!strcmp (op, "yramp")) {
    for (j = 0; j < Ny; j++) {
      for (i = 0; i < Nx; i++, out++, M1++) {
        *out = j;
      }
    }
  }
  
  if (V1[0].type == 'm') {
    free (V1[0].buffer[0].header.buffer);
    free (V1[0].buffer[0].matrix.buffer);
    free (V1[0].buffer);
  }

  clear_stack (V1);
  return (TRUE);

}

/*********************** fits copy header ***********************************/
int gfits_copy_matrix_info (Matrix *matrix1, Matrix *matrix2) {

  int i;

  /* copy all but the matrix */

  matrix2[0].unsign = matrix1[0].unsign;
  matrix2[0].bitpix = matrix1[0].bitpix;
  matrix2[0].size   = matrix1[0].size;
  matrix2[0].bzero  = matrix1[0].bzero;
  matrix2[0].bscale = matrix1[0].bscale;
  matrix2[0].Naxes  = matrix1[0].Naxes;
  for (i = 0; i < FT_MAX_NAXES; i++) 
    matrix2[0].Naxis[i] = matrix1[0].Naxis[i];


  return (TRUE);
}       
