/* file subsc.c  ana subscript routines */
#include <stdarg.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "ana_structures.h"
#include "defs.h"
#define	MIN(a,b) (((a)<(b))?(a):(b))
#define	MAX(a,b) (((a)>(b))?(a):(b))
 extern	struct sym_desc sym[];
 extern	int	temp_base;
 extern	char *find_sym_name(int iq);
 union	types_ptr { byte *b; short *w; int *l; float *f; double *d;};		
 extern	int	ana_type_size[];
 extern	int	ana_zero(), execute_error(), ana_byte(), ana_word();
 extern	int	ana_long(), ana_float(), ana_double();
 extern	int	string_scratch(), strarr_scratch();
 extern char	*strsave();
 static	int	result_sym;
 static	int	rear[10], loop[10], press[10], inc[10], lb[10];
 static	int	offset[10] ,comf[10];
/* 3/25/94  base needs to be an I*8, not an int, for 64 bit machines
 we assume that for 32 bit, long is I*4 and for 64 bit, it is I*8 */
 static long	base[10];	/* long is 64 bits on alpha and 64 bit Mac */
 static	int	scratch[10];
 static	int	ncomp, nbtriv, ind, nbdim, nd;
 /*----------------------------------------------------------------------- */
int ana_inserter(narg,ps)
 /* INSERT, array, sub, [i1,i2,...] */
 /* but currently limited to 1 and 2-D arrays */
 int narg,ps[];
 {
 int	iq,jq,nx,ny,nd1,nd2,type1,type2,ix,iy,nx2,ny2,nc,del,del2, nx2s;
 union	types_ptr p1, p2;
 struct	ahead	*h1, *h2;
 /* first must be an array */
 iq = ps[0];
 if ( sym[iq].class != 4 ) return execute_error(66);
 type1 = sym[iq].type;
 h1 = (struct ahead *) sym[iq].spec.array.ptr;
 nd1 = h1->ndim;
 if (nd1 > 2) return execute_error(110);
 nx = h1->dims[0];
 if ( nd1 > 1 ) ny = h1->dims[1]; else ny = 1;
 p1.l = (int *) ((char *)h1 + sizeof(struct ahead));
					 /* second must also be an array */
 iq = ps[1];
 if ( sym[iq].class != 4 ) return execute_error(66);
 type2 = sym[iq].type;
			 /* convert second to type of first if necessary */
 if (type1 != type2) { switch (type1) {
	 case 0:	iq = ana_byte(1,&iq);	break;
	 case 1:	iq = ana_word(1,&iq);	break;
	 case 2:	iq = ana_long(1,&iq);	break;
	 case 3:	iq = ana_float(1,&iq);	break;
	 case 4:	iq = ana_double(1,&iq);	break; } }
 h2 = (struct ahead *) sym[iq].spec.array.ptr;
 nd2 = h2->ndim;
 if (nd2 > 2) return execute_error(110);
 nx2 = nx2s = h2->dims[0];
 if ( nd2 > 1 ) ny2 = h2->dims[1]; else ny2 = 1;
 p2.l = (int *) ((char *)h2 + sizeof(struct ahead));
 if (narg >2) ix = int_arg( ps[2]); else ix = 0;
 if (narg >3) iy = int_arg( ps[3]); else iy = 0;
					 /* compute starting offset */
 /* 1/14/2001 - another bug when both x and y offset are negative */
 /* 7/17/2000 - serious bug, was not handling negative indices
 correctly and was clobbering arrays */
 /* 7/17/2000 - if ix or iy are negative, we now clip the second array
 if any part of it would fit in, also check for other out of bounds
 conditions, don't complain, just do nothing if completely out of
 bounds */
 if (ix >=0 ) {
  if (ix >= nx) return 1;
  iq = ix;
  jq = 0;
  del2 = 0;
  if (nx2 > (nx-ix) ) {
   /* extends over the right */
   del2 = nx2 - (nx - ix);
   nx2 = nx -ix;
  }
 } else {
  if (ix <= -nx2) return 1;
  jq = -ix;
  iq = 0;
  del2 = -ix;
  nx2 = nx2 + ix;
  if (nx2 > nx) {
   /* must extend over both sides */
   del2 = del2 + (nx2-nx);
   nx2 = nx;
  }
 }
 if (iy >= 0 ) {
  if (iy >= ny) return 1;
  /* no change to jq */
  iq = iq + iy*nx;
  if (ny2 > (ny-iy) ) ny2 = ny -iy;
 } else {
  if (iy <= -ny2) return 1;
  /* no change to iq */
  //jq = jq - iy * nx2;		/* bug (1/14/01) */
  jq = jq - iy * nx2s;		/* the fix (1/14/01) */
  ny2 = ny2 + iy;
  if (ny2 > ny) ny2 = ny;
 }
 del = nx - nx2;			/*leftovers per row */
 /* use only 4 switches for 4 sizes */
 /* printf("nc, ny2 = %d, %d\n", nc, ny2); */
 switch (type1) {
 case 0:
 p1.b += iq;
 p2.b += jq;
 while (ny2--) { nc = nx2; 
   while (nc--) *p1.b++ = *p2.b++;  p1.b += del; p2.b += del2; } break;
 case 1:
 p1.w += iq;
 p2.w += jq;
 while (ny2--) { nc = nx2; 
   while (nc--) *p1.w++ = *p2.w++;  p1.w += del; p2.w += del2; } break;
 case 2:
 case 3:
 p1.l += iq;
 p2.l += jq;
 while (ny2--) { nc = nx2; 
   while (nc--) *p1.l++ = *p2.l++;  p1.l += del; p2.l += del2; } break;
 case 4:
 p1.d += iq;
 p2.d += jq;
 while (ny2--) { nc = nx2; 
   while (nc--) *p1.d++ = *p2.d++;  p1.d += del; p2.d += del2; } break;
 }
 return 1;
 }
 /*------------------------------------------------------------------------- */
int ana_smap(narg,ps)
 /* convert type (and class) to a string without changing memory contents */
 int narg,ps[];
 {
 int n;
 union	types_ptr q1,q3;
 int	nsym, nd, j;
 struct	ahead	*h;
 nsym=ps[0];				/*only one argument */
 if (sym[nsym].class == 2 ) return nsym;	/*already a string */
 /*switch on the class */
 switch (sym[nsym].class)	{
 case 8:							/*scalar ptr*/
  n=1; q1.l = sym[nsym].spec.array.ptr;  break;
 case 1:								/*scalar */
  n=1; q1.l = &sym[nsym].spec.scalar.l;  break;
 case 4:								/*array */
  h = (struct ahead *) sym[nsym].spec.array.ptr;
  q1.l = (int *) ((char *)h + sizeof(struct ahead));
  nd = h->ndim;
  n = 1; for (j=0;j<nd;j++) n *= h->dims[j];	/*# of elements for nsym */
 break;
 default:	return execute_error(32);
 }
 n = n * ana_type_size[sym[nsym].type];
 result_sym = string_scratch( n );
 q3.l = (int *) sym[result_sym].spec.array.ptr;
 bcopy ( (char *) q1.b, (char *) q3.b, n );
 *(q3.b + n) = 0;
 return result_sym;
 }
 /*------------------------------------------------------------------------- */
int ana_bmap(narg,ps)
		/* convert type to byte without changing memory contents */
 int narg,ps[];
 {
 int	new;
 new = 0;
 return ana_gmap(narg,ps,new);
 }
 /*------------------------------------------------------------------------- */
int ana_wmap(narg,ps)
		/* convert type to byte without changing memory contents */
 int narg,ps[];
 {
 int	new;
 new = 1;
 return ana_gmap(narg,ps,new);
 }
 /*------------------------------------------------------------------------- */
int ana_lmap(narg,ps)
		/* convert type to long without changing memory contents */
 int narg,ps[];
 {
 int	new;
 new = 2;
 return ana_gmap(narg,ps,new);
 }
 /*------------------------------------------------------------------------- */
int ana_fmap(narg,ps)
		/* convert type to float without changing memory contents */
 int narg,ps[];
 {
 int	new;
 new = 3;
 return ana_gmap(narg,ps,new);
 }
 /*------------------------------------------------------------------------- */
int ana_dmap(narg,ps)
		/* convert type to double without changing memory contents */
 int narg,ps[];
 {
 int	new;
 new = 4;
 return ana_gmap(narg,ps,new);
 }
 /*------------------------------------------------------------------------- */
int ana_gmap(narg,ps,new)
		/* general part for map routines */
 int narg,ps[],new;
 {
 union	types_ptr q1,q3;
 int	nsym, nd, j, n, nn;
 int	type;
 struct	ahead	*h;
 nsym=ps[0];				/*only one argument */
 type = sym[nsym].type;
 /*check if already the desired type */
 if ( type == new ) return nsym;
 nd = 0;
 /*switch on the class */
 switch (sym[nsym].class)	{
 case 8:							/*scalar ptr*/
  n = ana_type_size[type]; q1.l = sym[nsym].spec.array.ptr;  break;
 case 1:								/*scalar */
  n = ana_type_size[type]; q1.l = &sym[nsym].spec.scalar.l;  break;
 case 2:							/*scalar ptr*/
  n = sym[nsym].spec.array.bstore-1; q1.l = sym[nsym].spec.array.ptr;  break;
 case 4:								/*array */
  h = (struct ahead *) sym[nsym].spec.array.ptr;
  q1.l = (int *) ((char *)h + sizeof(struct ahead));
  nd = h->ndim;
  n = h->dims[0] * ana_type_size[type];
 /* n must be a multiple of type size or the inner dimension is not
		 compatiable*/
  if ( (n % ana_type_size[new])  != 0 ) return execute_error(80);
  for (j=1;j<nd;j++) n *= h->dims[j];	/*# of elements for nsym */
  break;
 default:	return execute_error(32);
 }
		 /* n is # of bytes to transfer, get # of elements */
 nn = n/ ana_type_size[new];
 if ( nn <= 1 )  {	result_sym = scalar_scratch( new );
			 q3.l = &sym[result_sym].spec.scalar.l;
			 *q3.l = 0; }
		 /* more than 1 but not an array, make result an array */
 else { if ( nd == 0 ) { nd = 1;
	 result_sym = array_scratch( new, nd, &nn); }
		 /* input was an array, make duplicate and then change */
  else { result_sym = array_clone( nsym, type);
	 sym[result_sym].type = new; 
	 h = (struct ahead *) sym[result_sym].spec.array.ptr;
	 /*already know that inner dim is a proper multiple so safe to
		 do the following */
	 h->dims[0] = h->dims[0] * ana_type_size[type] / ana_type_size[new]; 
	 }
 h = (struct ahead *) sym[result_sym].spec.array.ptr;
 q3.l = (int *) ((char *)h + sizeof(struct ahead));
 }
 if (n>0) bcopy ( (char *) q1.b, (char *) q3.b, n );
 return result_sym;
 }
 /*------------------------------------------------------------------------- */
int ana_sub_arg(narg,ps)
 /* processing generalized subscripts */
 int narg,ps[];
 {
 struct	ahead	*h;
 int	i, *p, stat;
 i = 4;
 result_sym = array_scratch(2,1,&i);
 h = (struct ahead *) sym[result_sym].spec.array.ptr;
 p = (int *) ((char *)h + sizeof(struct ahead));
 sym[result_sym].class = SUBSC_DESC;		/*note class modification */
 if ( narg != 4 ) return execute_error(38);
 for (i=0;i<4;i++) {
   if (int_arg_stat(ps[i], p) != 1) return execute_error(39);  p++; }
 return result_sym;
 }
 /*------------------------------------------------------------------------- */
void rearr(x)
 /* rearrange subscript vectors */
 int	x[];
 {
 int	k;
 for (k=0;k<nd;k++) scratch[k] = x[ rear[k] ];
 for (k=0;k<nd;k++) x[k] = scratch[k];
 }
 /*------------------------------------------------------------------------- */
int ana_reverse(narg,ps)
						/* general array reverser */
 int narg,ps[];
 {
 union	types_ptr p1,p2;
 int	i, ns, j;
 int	nsym, iq, type, n, k, stride, mq;
 struct	ahead	*h;
 nsym = ps[0];				/*the target symbol is the first */
 type = sym[nsym].type;
				 /* only array and strings can be reversed */
 switch (sym[nsym].class)	{			/*switch on class */
 default: return execute_error(66);
 case STRING:							/*string case */
  p1.l = sym[nsym].spec.array.ptr;
  mq = sym[nsym].spec.array.bstore - 1;
  result_sym = string_scratch(mq);		/*for resultant string */
  p2.l = sym[result_sym].spec.array.ptr;
  simple_reverse(p1.l, p2.l, mq, 0);
  /* 4/20/95 bug found by L. Strous, no null at end of string */
  *(p2.b+mq) = 0;
  return	result_sym;
 case STRING_PTR:						/*strptr case */
  { char **p, **q, *q2;
  int	s_offset, dim[8];
  /* similar to the array case but we have to copy all the strings */
  h = (struct ahead *) sym[nsym].spec.array.ptr;
  p = (char **) ((char *) h + sizeof(struct ahead));
  nd = h->ndim;					/* number of dimensions */
  n = 1; for (j=0;j<nd;j++) { dim[j] = h->dims[j]; n *= h->dims[j]; }
  /* check if too many reverses specified */
  if (narg > (nd+1)) return execute_error(63);
  stride = 1;	/* one for pointer offsets */
  result_sym = strarr_scratch(nd, dim);
  h = (struct ahead *) sym[result_sym].spec.array.ptr;
  q = (char **) ((char *) h + sizeof(struct ahead));
  if (narg == 1)  {
  /* just back to front */
  p += n;
  while (n--) {
   q2 = *(--p);
   if (q2) *q = strsave(q2);
   q++;
  }
  
  return result_sym; }

  /* first set up as if there were no reversals */
  
  for (j=0;j<nd;j++) { inc[j] = stride;  stride *= h->dims[j];
	  loop[j] = lb[j] = h->dims[j]; }
  loop[nd] = lb[nd] = 0;	inc[nd] = stride;  base[0] = 0;
  /* loop over arguments, each must be an integer and less than nd */
  for (j=1;j<narg;j++) {
	  i = int_arg( ps[j] );
	  if (i < 0 || i >= nd ) return execute_error(81);
	  if (inc[i] < 0 ) return execute_error(82);
	  base[0] -= inc[i];
	  inc[i] = - inc[i];
	  iq = inc[i+1];
	  base[0] += abs(iq);
  }
  /* trickle up base */
  for (k=1;k<nd;k++) base[k] = base[0];
  j = inc[0];		/* inner loop increment */

  /* outer loops */
  //while (1) { s_offset = (byte *) base[1]; ns = lb[0];
  while (1) { s_offset = base[1]; ns = lb[0];

    /* inner loop */
    while (ns) {
	  q2 = *(p + s_offset);
	  if (q2) *q = strsave(q2);
	  q++;
 	  s_offset += j; ns--; }
    i = 1;
    /* process outer counts */
	    while (1) { loop[i] -= 1;
	    if ( loop[i] > 0 ) break;
	    if (loop[i] < 0) return result_sym;
	    i++; }
    base[i] += inc[i];
    while (i > 1) { i--; base[i] = base [i+1]; loop[i] = lb[i]; }
  }					/*end of outer loop while */
   /* note - returns are in the type switch cases above */
 }

 case ARRAY:							/*array case */
 h = (struct ahead *) sym[nsym].spec.array.ptr;
 p1.l = (int *) ((char *)h + sizeof(struct ahead));
 nd = h->ndim;					/* number of dimensions */
 n = 1; for (j=0;j<nd;j++) n *= h->dims[j]; 	/* # of elements for nsym */
 /* check if too many reverses specified */
 if (narg > (nd+1)) return execute_error(63);
 stride = ana_type_size[type];
 result_sym = array_clone( nsym ,type);		/* result same as input */
 h = (struct ahead *) sym[result_sym].spec.array.ptr;
 p2.l = (int *) ((char *)h + sizeof(struct ahead));
 if (narg == 1)  {	simple_reverse(p1.l, p2.l, n, type);
			 return result_sym; }
			 /* first set up as if there were no reversals*/
 
 for (j=0;j<nd;j++) { inc[j] = stride;  stride *= h->dims[j];
	 loop[j] = lb[j] = h->dims[j]; }
 loop[nd] = lb[nd] = 0;	inc[nd] = stride;  base[0] = (long) p1.l;
	 /* loop over arguments, each must be an integer and less than nd */
 for (j=1;j<narg;j++) {
	 i = int_arg( ps[j] );
	 if (i < 0 || i >= nd ) return execute_error(81);
	 if (inc[i] < 0 ) return execute_error(82);
	 base[0] -= inc[i];
	 inc[i] = - inc[i];
	 iq = inc[i+1];
	 base[0] += abs(iq);
 }
 /*trickle up base */
 for (k=1;k<nd;k++) base[k] = base[0];
 j = inc[0]/ana_type_size[type];		/*inner loop increment */
						 /* switch on source type */
 switch (type) {
 
 case 1:
 /*outer loops */
 while (1) { p1.w = (short *) base[1]; ns = lb[0];
 /* inner loop */
 while (ns) { *p2.w++ = *p1.w; p1.w += j; ns--; }
 i = 1;
 /*process outer counts */
	 while (1) { 
	 loop[i] -= 1;
	 if ( loop[i] > 0 ) break;
	 if (loop[i] < 0) return result_sym;
	 i++; }
 base[i] += inc[i];
 while (i > 1) { i--; base[i] = base [i+1]; loop[i] = lb[i]; }
 }					/*end of outer loop while */
 
 case 2:
 /*outer loops */
 while (1) { p1.l = (int *) base[1]; ns = lb[0];
 /* inner loop */
 while (ns) { *p2.l++ = *p1.l; p1.l += j; ns--; }
 i = 1;
 /*process outer counts */
	 while (1) { loop[i] -= 1;
	 if ( loop[i] > 0 ) break;
	 if (loop[i] < 0) return result_sym;
	 i++; }
 base[i] += inc[i];
 while (i > 1) { i--; base[i] = base [i+1]; loop[i] = lb[i]; }
 }					/*end of outer loop while */
 
 case 3:
 /*outer loops */
 while (1) { p1.f = (float *) base[1]; ns = lb[0];
 /* inner loop */
 while (ns) { *p2.f++ = *p1.f; p1.f += j; ns--; }
 i = 1;
 /*process outer counts */
	 while (1) { loop[i] -= 1;
	 if ( loop[i] > 0 ) break;
	 if (loop[i] < 0) return result_sym;
	 i++; }
 base[i] += inc[i];
 while (i > 1) { i--; base[i] = base [i+1]; loop[i] = lb[i]; }
 }					/*end of outer loop while */
 
 case 0:
 /*outer loops */
 while (1) { p1.b = (byte *) base[1]; ns = lb[0];
 /* inner loop */
 while (ns) { *p2.b++ = *p1.b; p1.b += j; ns--; }
 i = 1;
 /*process outer counts */
	 while (1) { loop[i] -= 1;
	 if ( loop[i] > 0 ) break;
	 if (loop[i] < 0) return result_sym;
	 i++; }
 base[i] += inc[i];
 while (i > 1) { i--; base[i] = base [i+1]; loop[i] = lb[i]; }
 }					/*end of outer loop while */
 
 case 4:
 /*outer loops */
 while (1) { p1.d = (double *) base[1]; ns = lb[0];
 /* inner loop */
 while (ns) { *p2.d++ = *p1.d; p1.d += j; ns--; }
 i = 1;
 /*process outer counts */
	 while (1) { loop[i] -= 1;
	 if ( loop[i] > 0 ) break;
	 if (loop[i] < 0) return result_sym;
	 i++; }
 base[i] += inc[i];
 while (i > 1) { i--; base[i] = base [i+1]; loop[i] = lb[i]; }
 }					/*end of outer loop while */
 }
 }						/* end of class switch */
 /*note - returns are in the type switch cases above */
 }
 /*------------------------------------------------------------------------- */
int simple_reverse(p1,p2,n,type)
 int	n, type;
 int	*p1, *p2;
 {
 union	types_ptr q1, q2;
 q1.l = p1;	q2.l = p2;
 switch (type) {
 case 0:	q1.b += n;  while (n--) *q2.b ++ = *(--q1.b); break;
 case 1:	q1.w += n;  while (n--) *q2.w ++ = *(--q1.w); break;
 case 2:	q1.l += n;  while (n--) *q2.l ++ = *(--q1.l); break;
 case 3:	q1.f += n;  while (n--) *q2.f ++ = *(--q1.f); break;
 case 4:	q1.d += n;  while (n--) *q2.d ++ = *(--q1.d); break;
 }
 return 1;
 }
 /*------------------------------------------------------------------------- */
int ana_subsc_func(narg,ps)
 /*decipher and extract subscripted values from a symbol, supports arrays,
 strings, associated variables, and function pointers */
 int narg,ps[];
 {
 union	types_ptr p1,p2,p3;
 int	i, ns, j;
 int	nsym, iq, type, n, stype, k, stride, len;
 byte	*ptr;
 struct	ahead	*h, *h2;
 nsym = ps[0];				/*the target symbol is the first */
 type = sym[nsym].type;
 switch (sym[nsym].class)	{			/*switch on class */
   default:
	   printf("symbol name: %s\n", find_sym_name(nsym) );
	   return execute_error(31);
   case ASSOC:						/*assoc. var. case */
 	  /* assoc variables come here but are handed off to files.c */
 	  return ana_assoc_input(narg, ps);
   case STRING_PTR:					/*strarr case */
	  /* process in a separate routine, return resulting string */
	  return strptr_sub(narg, ps);
   case SYM_ARR:						/*sym ptr case */
	  /* process in a separate routine, return result */
	  return symptr_sub(narg, ps);
   case 7:						/*func. ptr. case */
 	  /* actually, these don't exist, but may have to use someday */
 	  printf("function pointers not supported yet\n"); return 4;
   case STRING:						/*string case */
	  return string_sub(narg, ps);
   case ARRAY:						/*array case */
	  h = (struct ahead *) sym[nsym].spec.array.ptr;
 	  p1.l = (int *) ((char *)h + sizeof(struct ahead));
 	  nd = h->ndim;				/*number of dimensions */
 	  n = 1; for (j=0;j<nd;j++) n *= h->dims[j];
 	  break;
   case SYM_PTR:				/*symbol ptr (not sym array) case */
     { int sym_key, *ptable, nkeys;
     /* a pointer symbol, subscript indicates which of a list of symbols
     to return */
     /* this case only accepts a single subscript, do it all in line here */
     if (narg != 2) { return execute_error(63); }
     /* the second symbol must be a scalar which gets converted to int */
     if (int_arg_stat(ps[1], &sym_key) != 1) return -1;
     nkeys = sym[nsym].type;
     if (sym_key < 0 || sym_key >= nkeys) { return execute_error(35); }
     ptable = (int *) sym[nsym].spec.array.ptr;
     return	ptable[sym_key];
     }
     break;
  }				/*end of switch on class */
 /* actually, we get here only for class 4 (arrays) */
 if (narg <= 2 ) {			/*subscript is a single value */
   iq = ps[1];
   switch (sym[iq].class)	{
     case SUBSC_DESC: break;
     case SCAL_PTR: iq = class8_to_1(iq);
     case 1:	/*a scalar, simple, get a long (int) version */
     i = int_arg(iq);
     if ( i >= n || i < 0 ) {
 	   printf("array name: %s, index %d\n", find_sym_name(nsym), i );
   	   return execute_error(35); }
     return  create_sub_ptr(nsym,p1.l,i);

     case ARRAY:			/*an array subscript, only 1 allowed */
	    /*we create an array with structure of subscript array, the type
	    of main array, and containing elements of main array which
	    have subscripts equal to values of subscript array*/
     result_sym = array_clone(iq, type);
     h = (struct ahead *) sym[result_sym].spec.array.ptr;
     p3.l = (int *) ((char *)h + sizeof(struct ahead));
     nd = h->ndim;
     ns = 1; for (j=0;j<nd;j++) ns *= h->dims[j]; 	/*# of elements for result */
     h = (struct ahead *) sym[iq].spec.array.ptr;
     p2.l = (int *) ((char *)h + sizeof(struct ahead));
     stype = sym[iq].type;
     while (ns) {					/*loop over result range */
     /*get long version of index */
     /*might be faster (but more memory intensive) to make a temp copy of long
	    type, a future experiment */
     switch (stype) {
     case 2:	i = *p2.l++; break;
     case 0:	i = (int) *p2.b++; break;
     case 1:	i = (int) *p2.w++; break;
     case 3:	i = (int) *p2.f++; break;
     case 4:	i = (int) *p2.d++; break;
     }

    /* 10/11/92 */
    /* to make lookup tables work better, force out of range values to 0 or n-1 */
    /* if ( i >= n || i < 0 ) return execute_error(35);*/
    if ( i < 0 ) i = 0;	if ( i >= n ) i = n-1;
    /*apply and load into result */	
     switch (type) {
     case 0:	*p3.b++ = (byte) *( (p1.b) + i ); break;	
     case 1:	*p3.w++ = (short) *( (p1.w) + i ); break;	
     case 2:	*p3.l++ = (int) *( (p1.l) + i ); break;	
     case 3:	*p3.f++ = (float) *( (p1.f) + i ); break;	
     case 4:	*p3.d++ = (double) *( (p1.d) + i ); break;
     }	
     ns--; }
     return result_sym;
   }
 }
							 /*general case */
	 /*several subscripts, the count (narg-1) can't be more than nd */
 if (narg > (nd+1) ) return execute_error(36);
 stride = ana_type_size[type];
 ptr = p1.b;
 /* printf("ptr, p1.f, *p1.f %d %d %f\n", ptr, p1.f, *p1.f);*/
 for (k = 1; k < narg; k++) {			/*loop over arguments */
 iq = ps[k];	j = h->dims[k-1];
 switch (sym[iq].class)	{	/*switch on class of this subscript */
 case SCAL_PTR: iq = class8_to_1(iq);
 case 1: /*scalar*/
  i = int_arg(iq);
  if ( i >= j || i < 0 ) {
 	printf("array name: %s, index %d\n", find_sym_name(nsym), i );
   	return execute_error(35); }
  ptr = ptr + (i*stride);	stride = stride*j;	break;
 case SUBSC_DESC: /*special subscript type (only allowed type here), we do not
	 return to the k loop, instead a new loop is started where the
	 old one left off, a return is executed before we go back */
  nd = nd - k + 1;		/* decrease nd by amount already done */
  for (j=0;j<10;j++) { loop[j] = 0; rear[j] = -1; }		/*set ups */
  loop[0] = 0;

  base[0] = (long) ptr;
  /* printf("ptr, base[0] = %d %d\n", ptr, base[0]);*/
  /* printf(" *(float *)ptr = %f\n",  *(float *)ptr ); */
  ncomp=nbtriv=nbdim=ind=0;
  /*ind is dimension count in non-simple mode */
  /*nbtriv is the # of trivial dimensions */
  while (k < narg) {	/*loop over remaining dimensions */
  press[ind] = 0;
  inc[ind] = stride;	/*save stride in inc array for each dim. */
  j = h->dims[k-1];
  stride = stride*j;
  iq = ps[k];		/*current index, check if simple */
  if (sym[iq].class == 8) iq = class8_to_1(iq);
  if (sym[iq].class == 1) {  i = int_arg(iq);		/* simple case */
	  /*a simple one in the list */
	  if ( i >= j || i < 0 ) {
 		printf("array name: %s, index %d\n", find_sym_name(nsym), i );
   		return execute_error(35); }
	  loop[ind] = lb[ind] = 1;	comf[ind] = 0;
	  offset[ind] = i*inc[ind]; base[0] += offset[ind];
	  if ( rear[ind] >= 0 ) return execute_error(37);
	  rear[ind] = ind;	nbtriv += 1;}
 else {			/*non-simple case, check if the allowed type */
 if (sym[iq].class != SUBSC_DESC) return execute_error(39);
 /*the class SUBSC_DESC symbol contains 4 integers, pull them out */
 h2 = (struct ahead *) sym[iq].spec.array.ptr;
 p2.l = (int *) ((char *)h2 + sizeof(struct ahead));
 i = *p2.l++;	/*first one is the starting index */ 
 if ( i >= j || i < 0 ) {
 	printf("array name: %s, index %d\n", find_sym_name(nsym), i );
   	return execute_error(35); }
 len = *p2.l++;	/*second is the ending index */
 /* printf("start = %d, len = %d\n", i, len); */
 if ( len < 0 ) len = j - i; else  {
	 if ( len >= j || len < i ) {
 		printf("array name: %s, index %d\n", find_sym_name(nsym), len );
   		return execute_error(35); }
	 len = len - i + 1; }
 loop[ind] = lb[ind] = len;
 if (comf[ind] = *p2.l++) ncomp++;
 offset[ind] = i*inc[ind]; base[0] += offset[ind];
 /*now check the rearrange (fourth quantity) */
 if ( len == 1 || comf[ind] != 0 ) {	/*must be trivial, check it */
	  if ( rear[ind] >= 0 ) return execute_error(37);
	  rear[ind] = ind;	nbtriv += 1;}
 else {	len = *p2.l++;	/*re-use len variable for rearrange target */
	 if ( len < 0 ) len = nbdim;
	 len = len + nbtriv; if ( rear[len] >= 0 ) return execute_error(37);
	 rear[len] = ind;	nbdim += 1; }
 }
 ind += 1;
 k++; }
 /*still in case SUBSC_DESC section, looping over subscripts finished */
 if ( ind != nd ) {
 printf("dimension count error, outer dimensions will be assumed 0\n");
 printf("nd, ind = %d, %d\n",nd, ind); nd = ind; }
 for (k=0;k<nd;k++) if (rear[k] < 0 ) return execute_error(37);
 /*trickle up base */
 for (k=1;k<=nd;k++) base[k] = base[0];
 /*rearrange everybody */
 (void) rearr(inc);	(void) rearr(loop);	(void) rearr(lb);	
 /*setup output array, get actual dimensions in scratch */
 j = 0;
 for (k=0;k<nd;k++) if (loop[k] > 1 && comf[k] == 0 ) scratch[j++]=loop[k];
 j = inc[0]/ana_type_size[type];		/*inner loop increment */
 /*the result will always be float if there is compression, otherwise
 same type as source */
 if (ncomp) {
 /*might have compressed to only one element, if so, make a scalar */
	 if (nbdim) {			/*array case */
	 result_sym = array_scratch( MAX(3,type), nbdim, scratch);
				 /*call zero to zero this array */
	 if ( ana_zero(1,&result_sym) != 1 ) return execute_error(3);
	 } else {			/*scalar result */
	 result_sym = scalar_scratch( MAX(3,type));
	 sym[result_sym].spec.scalar.d = 0.0;
	 p2.d = &sym[result_sym].spec.scalar.d; }
 } else {
  /* non-compression, if reduced to 1 element, make a scalar here also */
  ns = 1; for (k=0;k<nbdim;k++) ns *= scratch[k];
  if (ns == 1) { result_sym = scalar_scratch(type);
	 sym[result_sym].spec.scalar.d = 0.0;
	 p2.d = &sym[result_sym].spec.scalar.d; }
		else result_sym = array_scratch( type, nbdim, scratch);
  }
 if (sym[result_sym].class == 4) {
	 h = (struct ahead *) sym[result_sym].spec.array.ptr;
	 p2.l = (int *) ((char *)h + sizeof(struct ahead)); }
 /*branch on whether there is compression */
 if (ncomp) {					/*compression case */
 /*set up compression backlashes */
 press[0] = 1;
 for (k=1;k<nd;k++) {
	 /*look at previous one, was it compressed ? */
	 if ( comf[ rear[k-1] ]) press[k] = press[k-1];	/* yes */
	 else { press[k] = lb[k-1] * press[k-1]; press[k-1]=0; } /* no */
 }	/*still need to check the last one */
 if (comf[nd-1] == 0) press[nd-1] = 0;
 /* switch on source type */
 switch (type) {
 
 case 1:
 /*outer loops */
 while (1) { p1.w = (short *) base[1]; ns = lb[0]; i = 1 - press[0];
 /* inner loop */
 while (ns) { *p2.f += (float) *p1.w; p1.w += j; p2.f += i; ns--; }
 p2.f += press[0];	i = 1;
 /*process outer counts */
	 while (1) { loop[i] -= 1;
	 if ( loop[i] > 0 ) break;
	 if (loop[i] < 0) return result_sym;
	 i++; }
 base[i] += inc[i]; p2.f -= press[i];
 while (i > 1) { i--; base[i] = base [i+1]; loop[i] = lb[i]; }
 }					/*end of outer loop while */
 case 2:
 /*outer loops */
 while (1) { p1.l = (int *) base[1]; ns = lb[0]; i = 1 - press[0];
 /* inner loop */
 while (ns) { *p2.f += (float) *p1.l; p1.l += j; p2.f += i; ns--; }
 p2.f += press[0];	i = 1;
 /*process outer counts */
	 while (1) { loop[i] -= 1;
	 if ( loop[i] > 0 ) break;
	 if (loop[i] < 0) return result_sym;
	 i++; }
 base[i] += inc[i]; p2.f -= press[i];
 while (i > 1) { i--; base[i] = base [i+1]; loop[i] = lb[i]; }
 }					/*end of outer loop while */
 
 case 3:
 /*outer loops */
 /*printf("p2.f = %d\n",p2.f );*/
 while (1) { p1.f = (float *) base[1]; ns = lb[0]; i = 1 - press[0];
 /* inner loop */
 while (ns) { 
 /*printf("in loop, p2.f = %d, p1.f = %d\n",p2.f, p1.f );*/
 *p2.f += (float) *p1.f; p1.f += j; p2.f += i; ns--; }
 p2.f += press[0];	i = 1;
 /*process outer counts */
	 while (1) { loop[i] -= 1;
	 if ( loop[i] > 0 ) break;
	 if (loop[i] < 0) return result_sym;
	 i++; }
 base[i] += inc[i]; p2.f -= press[i];
 while (i > 1) { i--; base[i] = base [i+1]; loop[i] = lb[i]; }
 }					/*end of outer loop while */
 
 case 0:
 /*outer loops */
 while (1) { p1.b = (byte *) base[1]; ns = lb[0]; i = 1 - press[0];
 /* inner loop */
 while (ns) { *p2.f += (float) *p1.b; p1.b += j; p2.f += i; ns--; }
 p2.f += press[0];	i = 1;
 /*process outer counts */
	 while (1) { loop[i] -= 1;
	 if ( loop[i] > 0 ) break;
	 if (loop[i] < 0) return result_sym;
	 i++; }
 base[i] += inc[i]; p2.f -= press[i];
 while (i > 1) { i--; base[i] = base [i+1]; loop[i] = lb[i]; }
 }					/*end of outer loop while */
 
 case 4:
 /*outer loops */
 while (1) { p1.d = (double *) base[1]; ns = lb[0]; i = 1 - press[0];
 /* inner loop */
 while (ns) { *p2.d += *p1.d; p1.d += j; p2.d += i; ns--; }
 p2.d += press[0];	i = 1;
 /*process outer counts */
	 while (1) { loop[i] -= 1;
	 if ( loop[i] > 0 ) break;
	 if (loop[i] < 0) return result_sym;
	 i++; }
 base[i] += inc[i]; p2.d -= press[i];
 while (i > 1) { i--; base[i] = base [i+1]; loop[i] = lb[i]; }
 }					/*end of outer loop while */
 }	/*end of switch in compression case */
 }					/*end of compression case */
 
 
 else	{				/*no compression case */
 /* switch on source type */
 switch (type) {
 
 case 1:
 /*outer loops */
 while (1) { p1.w = (short *) base[1]; ns = lb[0];
 /* inner loop */
 while (ns) { *p2.w++ = *p1.w; p1.w += j; ns--; }
 i = 1;
 /*process outer counts */
	 while (1) { 
	 loop[i] -= 1;
	 if ( loop[i] > 0 ) break;
	 if (loop[i] < 0) return result_sym;
	 i++; }
 base[i] += inc[i];
 while (i > 1) { i--; base[i] = base [i+1]; loop[i] = lb[i]; }
 }					/*end of outer loop while */
 
 case 2:
 /*outer loops */
 while (1) { p1.l = (int *) base[1]; ns = lb[0];
 /* inner loop */
 while (ns) { *p2.l++ = *p1.l; p1.l += j; ns--; }
 i = 1;
 /*process outer counts */
	 while (1) { loop[i] -= 1;
	 if ( loop[i] > 0 ) break;
	 if (loop[i] < 0) return result_sym;
	 i++; }
 base[i] += inc[i];
 while (i > 1) { i--; base[i] = base [i+1]; loop[i] = lb[i]; }
 }					/*end of outer loop while */
 
 case 3:
 /*outer loops */
 /* printf("p2.f = %d\n",p2.f ); */
 while (1) { p1.f = (float *) base[1]; ns = lb[0];
 /* inner loop */
 while (ns) { 
 /* printf("noncompress, in loop, p2.f = %d, p1.f = %d\n",p2.f, p1.f ); */
 /* printf("*p2.f =  %f\n", *p2.f); */
 /* printf("*p1.f =  %f\n", *p1.f); */
 *p2.f++ = *p1.f; p1.f += j; ns--; }

 i = 1;
 /*process outer counts */
	 while (1) { loop[i] -= 1;
	 if ( loop[i] > 0 ) break;
	 if (loop[i] < 0) return result_sym;
	 i++; }
 base[i] += inc[i];
 while (i > 1) { i--; base[i] = base [i+1]; loop[i] = lb[i]; }
 }					/*end of outer loop while */
 
 case 0:
 /*outer loops */
 while (1) { p1.b = (byte *) base[1]; ns = lb[0];
 /* inner loop */
 while (ns) { *p2.b++ = *p1.b; p1.b += j; ns--; }
 i = 1;
 /*process outer counts */
	 while (1) { loop[i] -= 1;
	 if ( loop[i] > 0 ) break;
	 if (loop[i] < 0) return result_sym;
	 i++; }
 base[i] += inc[i];
 while (i > 1) { i--; base[i] = base [i+1]; loop[i] = lb[i]; }
 }					/*end of outer loop while */
 
 case 4:
 /*outer loops */
 while (1) { p1.d = (double *) base[1]; ns = lb[0];
 /* inner loop */
 while (ns) { *p2.d++ = *p1.d; p1.d += j; ns--; }
 i = 1;
 /*process outer counts */
	 while (1) { loop[i] -= 1;
	 if ( loop[i] > 0 ) break;
	 if (loop[i] < 0) return result_sym;
	 i++; }
 base[i] += inc[i];
 while (i > 1) { i--; base[i] = base [i+1]; loop[i] = lb[i]; }
 }					/*end of outer loop while */
 }				/*end of switch in compression case */
 }				/*end of if (ncomp) condition */
 }				/*end of switch */
 }	/*we get here only for a series of scalar subscripts, point to a
	 single value */
 return  create_sub_ptr(nsym,ptr,0);
 }						/* end of ana_subsc */
 /*------------------------------------------------------------------------- */
int string_sub(narg,ps)	/* subscripts for strings, called by ana_subc */
 int narg,ps[];
 {
 /* this is a subset of cases for arrays since strings are 1-D */
 int	iq, n, result_sym, j, nd, ns, nsym, i;
 char	*p, *q;
 struct	ahead	*h;
 union	types_ptr p2;
 if ( narg != 2 ) return execute_error(87);
 nsym = ps[0];
 p = (char *) sym[nsym].spec.array.ptr;
 n = sym[nsym].spec.array.bstore - 1;
 iq = ps[1];					/* only one subscript */
  switch (sym[iq].class)	{
  case SCAL_PTR: iq = class8_to_1(iq);
  case SCALAR:	/*a scalar, simple, get a long (int) version */
  i = int_arg(iq);
  if ( i >= n || i < 0 ) {
 	printf("string name: %s, index %d\n", find_sym_name(nsym), i );
   	return execute_error(35); }
 /* pointers don't work too well for strings because various parts of code
	 expect strings to be null terminated, so for now we just create
	 a single element string, this means that an element can't be loaded
	 in a read statement like an array element */
 /* return  create_str_sub_ptr(nsym,p,i); */
  result_sym = string_scratch(1);
  q = (char *) sym[result_sym].spec.array.ptr;
  *q++ = *(p + i);	*q = 0;
  return result_sym;

  case ARRAY:
  /* a string is created with length = # elements in array */
  iq = ana_long(1, &iq);				/* get a long version */
  h = (struct ahead *) sym[iq].spec.array.ptr;
  p2.l = (int *) ((char *)h + sizeof(struct ahead));
  nd = h->ndim;
  ns = 1; for (j=0;j<nd;j++) ns *= h->dims[j]; 	/*# of elements for result */
  result_sym = string_scratch(ns);
  q = (char *) sym[result_sym].spec.array.ptr;
  while (ns) {
  i = *p2.l++;
  if ( i >= n || i < 0 ) {
 	printf("string name: %s, index %d\n", find_sym_name(nsym), i );
   	return execute_error(35); }
  *q++ = *( (p++) + i );
  ns--; }	*q = 0;
  return result_sym;

  case SUBSC_DESC:		/* not all aspects supported for strings */
	/*the class SUBSC_DESC symbol contains 4 integers, pull them out */
  h = (struct ahead *) sym[iq].spec.array.ptr;
  p2.l = (int *) ((char *)h + sizeof(struct ahead));
  i = *p2.l++;	/*first one is the starting index */ 
  if ( i >= n || i < 0 ) {
 	printf("string name: %s, index %d\n", find_sym_name(nsym), i );
   	return execute_error(35); }
  ns = *p2.l++;	/*second is the ending index */
  if ( ns < 0 ) ns = n - i; else  {
	 if ( ns >= n || ns < i ) {
 		printf("string name: %s, index %d\n", find_sym_name(nsym), ns );
   		return execute_error(35); }
	 ns = ns - i + 1; }
		 /* if the compression or rearrange are set, flag an error */
 
  if (*p2.l || (*(p2.l + 1)) >= 0 ) return execute_error(88);
  result_sym = string_scratch(ns);
  p += i;
  q = (char *) sym[result_sym].spec.array.ptr;
  while (ns) {
  *q++ = *p++;
  ns--; }	*q = 0;
  return result_sym;
  }
 }
 /*------------------------------------------------------------------------- */
int strarr_error(i, nsym)
 int	i, nsym;
 {
  printf("string array name: %s, index %d\n", find_sym_name(nsym), i );
  return execute_error(35);
 }
 /*------------------------------------------------------------------------- */
int strptr_sub(narg,ps)		/* extract a string from a string array */
 int narg,ps[];
 {
 /* extracts string(s) from the string array given the index,
 makes a copy and puts into a temporary symbol */
 int	iq, result_sym, j, nd, ns, nsym, nelem, s_offset, dim[10], *pint;
 int	i, n, k, stride, *pd;
 char	*p1, *p2, **p, **q, *q2, *ptr;
 struct	ahead	*h;
 /* if ( narg != 2 ) return execute_error(127); */
 /* upgrade to more than 1-D 9/22/97 */
 nsym = ps[0];
 h = (struct ahead *) sym[nsym].spec.array.ptr;
 p = (char **) ((char *)h + sizeof(struct ahead));
 nd = h->ndim;
 nelem = 1; for (j=0;j<nd;j++) nelem *= h->dims[j];	/*# of elements */

 if (narg <= 2 ) {			/*subscript is a single value */
 iq = ps[1];					/* only one subscript */
  switch (sym[iq].class)	{
  case SCAL_PTR: iq = class8_to_1(iq);
  case SCALAR:	/*a scalar, simple, get a long (int) version */
   if (int_arg_stat(ps[1], &s_offset) != 1) return execute_error(127);
   /* check if offset in range */
   if (s_offset < 0 || s_offset >= nelem)  return execute_error(127);
   /* get the pointer location */
   p =  p + s_offset;
   p2 = *p;
   if (p2) {
   ns = strlen(p2);
   result_sym = string_scratch(ns);
   strcpy( (char *) sym[result_sym].spec.array.ptr, p2);
   } else {
   /* a null pointer to the string, generate a string of 0 length */
   result_sym = string_scratch(0);
   * ( (char *) sym[result_sym].spec.array.ptr) = 0;
   }
   return result_sym;

 case ARRAY:
  /* array subscript, create a new string array with selected elements */
  iq = ana_long(1, &iq);			/* get a long version */
  h = (struct ahead *) sym[iq].spec.array.ptr;
  pint = (int *) ((char *)h + sizeof(struct ahead));
  nd = h->ndim;
  ns = 1; for (j=0;j<nd;j++) ns *= h->dims[j]; 	/*# of elements for result */
  /* get a string array of this size */
  result_sym = strarr_scratch(h->ndim, h->dims);
  q = (char **)((char *)sym[result_sym].spec.array.ptr+sizeof(struct ahead));
  while (ns--) {
  i = *pint++;	/* this is the desired index to extract */
  if ( i >= nelem || i < 0 ) {
 	printf("string array name: %s, index %d\n", find_sym_name(nsym), i );
   	return execute_error(127); }
  q2 = *(p+i);
  if (q2) *q = strsave(q2);
  q++;
  }
  return result_sym;
  case SUBSC_DESC:	/* not all aspects supported for string arrays */
	/*the class SUBSC_DESC symbol contains 4 integers, pull them out */
  h = (struct ahead *) sym[iq].spec.array.ptr;
  pint = (int *) ((char *)h + sizeof(struct ahead));
  i = *pint++;	/*first one is the starting index */ 
  ns = *pint++;	/*second is the ending index */
  if ( i >= nelem || i < 0 ) {
 	printf("string array name: %s, index %d\n", find_sym_name(nsym), i );
   	return execute_error(127); }
  if ( ns < 0 ) ns = nelem - i; else  {
	 if ( ns >= nelem || ns < i ) {
 	printf("string array name: %s, index %d\n", find_sym_name(nsym), i );
   	return execute_error(127); }
	 ns = ns - i + 1; }
	/* if the compression or rearrange are set, flag an error */
 
  if (*pint || (*(pint + 1)) >= 0 ) return execute_error(88);

  /* get a string array of this size */
  dim[0]=ns;	result_sym = strarr_scratch(1, dim);

  q = (char **)((char *)sym[result_sym].spec.array.ptr+sizeof(struct ahead));
  p = p + i;
  while (ns--) {
  q2 = *p++;
  if (q2) *q = strsave(q2);
  q++;
  }
  return result_sym;
  }
  } else {

 /* must be multi-D */
 /*several subscripts, the count (narg-1) can't be more than nd */
 if (narg > (nd+1) ) return execute_error(36);

 stride = 1;
 s_offset = 0;
 /* here we are manipulating pointers to strings, so use pointers rather than
 the addresses used for real arrays, hence stride starts as 1
 we finally get the offsets for the string pointers */
 for (k = 1; k < narg; k++) {			/*loop over arguments */
 iq = ps[k];	j = h->dims[k-1];
 switch (sym[iq].class)	{	/*switch on class of this subscript */
 case SCAL_PTR: iq = class8_to_1(iq);
 case SCALAR: 				/*scalar*/
  i = int_arg(iq);
  if ( i >= j || i < 0 ) return strarr_error(i, nsym);
  s_offset = s_offset + (i*stride);	stride = stride*j;
  /* printf("present s_offset = %d\n", s_offset); */
  break;
 case SUBSC_DESC: /*special subscript type (only allowed type here), we do not
	 return to the k loop, instead a new loop is started where the
	 old one left off, a return is executed before we go back */
  nd = nd - k + 1;		/* decrease nd by amount already done */
  /* this is similar to the array case but we are dealing with string pointers
  and compression is not supported (not sensible) */

  for (j=0;j<10;j++) { loop[j] = 0; rear[j] = -1; }		/*set ups */
  /* loop[0] = 0; */

  /* base[0] = (long) ptr; */
  base[0] = s_offset;	/* here as offsets from start, not absolute pointers */
  /* printf("ptr, base[0] = %d %d\n", ptr, base[0]);*/
  /* printf(" *(float *)ptr = %f\n",  *(float *)ptr ); */
  ncomp = nbtriv = nbdim = ind = 0;
  /* ind is dimension count in non-simple mode */
  /* nbtriv is the # of trivial dimensions */
  while (k < narg) {	/* loop over remaining dimensions */
  press[ind] = 0;
  inc[ind] = stride;	/* save stride in inc array for each dim. */
  j = h->dims[k-1];
  stride = stride*j;
  iq = ps[k];		/* current index, check if simple */
  if (sym[iq].class == 8) iq = class8_to_1(iq);
  if (sym[iq].class == 1) {  i = int_arg(iq);		/* simple case */
	  /* a simple one in the list */
	  if ( i >= j || i < 0 ) return strarr_error(i, nsym);
	  loop[ind] = lb[ind] = 1;	comf[ind] = 0;
	  offset[ind] = i*inc[ind]; base[0] += offset[ind];
	  if ( rear[ind] >= 0 ) return execute_error(37);
	  rear[ind] = ind;	nbtriv += 1;}
 else {			/*non-simple case, check if the allowed type */
 int	len;
 if (sym[iq].class != SUBSC_DESC) return execute_error(39);
 /* the class SUBSC_DESC symbol contains 4 integers, pull them out */
 pd = (int *) ((char *)sym[iq].spec.array.ptr + sizeof(struct ahead));
 i = *pd++;	/* first one is the starting index */ 
 if ( i >= j || i < 0 ) return strarr_error(i, nsym);
 len = *pd++;	/* second is the ending index */
 /* printf("start = %d, len = %d\n", i, len); */
 if ( len < 0 ) len = j - i; else  {
   if ( len >= j || len < i ) return strarr_error(i, nsym);
	 len = len - i + 1; }
 loop[ind] = lb[ind] = len;
 /* keep the compression info here in case we figure out a sensible thing
 to do with it */
 if (comf[ind] = *pd++) ncomp++;
 offset[ind] = i*inc[ind]; base[0] += offset[ind];
 /*now check the rearrange (fourth quantity) */
 if ( len == 1 || comf[ind] != 0 ) {	/*must be trivial, check it */
	  if ( rear[ind] >= 0 ) return execute_error(37);
	  rear[ind] = ind;	nbtriv += 1;}
 else {	len = *pd++;	/*re-use len variable for rearrange target */
	 if ( len < 0 ) len = nbdim;
	 len = len + nbtriv; if ( rear[len] >= 0 ) return execute_error(37);
	 rear[len] = ind;	nbdim += 1; }
 }
 ind += 1;
 k++; } 
 /* still in case SUBSC_DESC section, looping over subscripts finished */
 if ( ind != nd ) {
 printf("dimension count error, outer dimensions will be assumed 0\n");
 printf("nd, ind = %d, %d\n",nd, ind); nd = ind; }
 for (k=0;k<nd;k++) if (rear[k] < 0 ) return execute_error(37);
 /*trickle up base */
 for (k=1;k<=nd;k++) base[k] = base[0];
 /*rearrange everybody */
 (void) rearr(inc);	(void) rearr(loop);	(void) rearr(lb);	
 /*setup output array, get actual dimensions in scratch */
 j = 0;
 for (k=0;k<nd;k++) if (loop[k] > 1 && comf[k] == 0 ) scratch[j++]=loop[k];
 j = inc[0];		/*inner loop increment */
 if (ncomp) {
 /* compression (summation) not supported here */
	return execute_error(128); 
 } else {
  /* non-compression (always at present), if reduced to a single string,
  we just return the string rather than a string array */
  ns = 1; for (k=0;k<nbdim;k++) ns *= scratch[k];
  if (ns == 1) {
  	/* just get offset and break to single string stuff at end */
	s_offset = base[0];  break;
  } else {
  result_sym = strarr_scratch(nbdim, scratch);
  q = (char **)((char *)sym[result_sym].spec.array.ptr+sizeof(struct ahead));
  }
 
 while (1) { s_offset = base[1]; ns = lb[0];
 /* inner loop */
 while (ns) {
	q2 = *(p + s_offset);
	if (q2) *q = strsave(q2);
	q++;
 	s_offset += j; ns--; }
 i = 1;
 /*process outer counts */
	 while (1) { loop[i] -= 1;
	 if ( loop[i] > 0 ) break;
	 if (loop[i] < 0) return result_sym;
	 i++; }
 base[i] += inc[i];
 while (i > 1) { i--; base[i] = base [i+1]; loop[i] = lb[i]; }
 }					/*end of outer loop while */
 }
 }				/*end of switch */
 }
 /* we get here only for a series of scalar subscripts or a single string */
 /* check if offset in range */
 if (s_offset < 0 || s_offset >= nelem)  return execute_error(127);
 /* get the pointer location */
 /* printf("s_offset = %d\n", s_offset); */
 p =  p + s_offset;
 p2 = *p;
 if (p2) {
 ns = strlen(p2);
 result_sym = string_scratch(ns);
 strcpy( (char *) sym[result_sym].spec.array.ptr, p2);
 } else {
 /* a null pointer to the string, generate a string of 0 length */
 result_sym = string_scratch(0);
 * ( (char *) sym[result_sym].spec.array.ptr) = 0;
 }
 return result_sym;
 }
 }
 /*------------------------------------------------------------------------- */
int symptr_sub(narg,ps)		/* extract a symbol from a symbol array */
 int narg,ps[];
 {
 /* extracts symbol(s) from the symbol array given the index,
 makes a copy and puts into a temporary symbol or symbol array.
 Generally this type of action should be avoided since it means duplicating
 the symbols and their values. Routines that use a component of the symbol
 array should be used instead to avoid the copying (when possible). */
 int	iq, result_sym, j, nd, ns, nsym, nelem, s_offset, dim[10], *pint;
 int	i, n, k, stride, *pd, stat;
 //char	*p1, *p2, **p, **q, *q2, *ptr;
 struct sym_desc *p, *q;
 struct	ahead	*h;
 //printf("symptr_sub\n");
 nsym = ps[0];
 h = (struct ahead *) sym[nsym].spec.array.ptr;
 p = (struct sym_desc *) ((char *)h + sizeof(struct ahead));
 nd = h->ndim;
 nelem = 1; for (j=0;j<nd;j++) nelem *= h->dims[j];	/*# of elements */

 if (narg <= 2 ) {			/*subscript is a single value */
 iq = ps[1];				/* only one subscript */
  switch (sym[iq].class)	{	/* switch on class of subscript */
  case SCAL_PTR: iq = class8_to_1(iq);
  case SCALAR:	/* a scalar, simple, get a long (int) version */
   if (int_arg_stat(ps[1], &s_offset) != 1) return execute_error(127);
   /* check if offset in range */
   if (s_offset < 0 || s_offset >= nelem)  return execute_error(127);
   /* get the pointer location */
   p =  p + s_offset;
   /* get a scratch symbol, note that the result cannot be a symbol array */
   result_sym = find_next_temp();
   stat = symbol_clone_via_ptr(&sym[result_sym], p, 1); /* note mode = 1 here */
   if (stat != 1) return stat;
   return result_sym;

  case ARRAY:
   /* array subscript, create a new symbol array with selected elements */
   iq = ana_long(1, &iq);			/* get a long version */
   h = (struct ahead *) sym[iq].spec.array.ptr;
   pint = (int *) ((char *)h + sizeof(struct ahead));
   nd = h->ndim;
   ns = 1; for (j=0;j<nd;j++) ns *= h->dims[j]; 	/*# of elements for result */
   /* get a symbol array of this size */
   result_sym = symarr_scratch(h->ndim, h->dims);
   q = (struct sym_desc *)((char *)sym[result_sym].spec.array.ptr+sizeof(struct ahead));
   while (ns--) {
   i = *pint++;	/* this is the desired index to extract */
   if ( i >= nelem || i < 0 ) {
 	 printf("symbol array name: %s, index %d\n", find_sym_name(nsym), i );
   	 return execute_error(127); }
    stat = symbol_clone_via_ptr(q, p+i, 1);
    if (stat != 1) return stat;
    q++;
   }
   return result_sym;

  case SUBSC_DESC:	/* not all aspects supported for symbol arrays */
	/*the class SUBSC_DESC symbol contains 4 integers, pull them out */
   h = (struct ahead *) sym[iq].spec.array.ptr;
   pint = (int *) ((char *)h + sizeof(struct ahead));
   i = *pint++;	/*first one is the starting index */ 
   ns = *pint++;	/*second is the ending index */
   if ( i >= nelem || i < 0 ) {
 	 printf("string array name: %s, index %d\n", find_sym_name(nsym), i );
   	 return execute_error(127); }
   if ( ns < 0 ) ns = nelem - i; else  {
	  if ( ns >= nelem || ns < i ) {
 	 printf("string array name: %s, index %d\n", find_sym_name(nsym), i );
   	 return execute_error(127); }
	  ns = ns - i + 1; }
	 /* if the compression or rearrange are set, flag an error */

   if (*pint || (*(pint + 1)) >= 0 ) return execute_error(88);

   /* get a string array of this size */
   dim[0]=ns;	result_sym = symarr_scratch(1, dim);

   q = (struct sym_desc *)((char *)sym[result_sym].spec.array.ptr+sizeof(struct ahead));
   p = p + i;
   while (ns--) {
   stat = symbol_clone_via_ptr(q, p, 1);
   if (stat != 1) return stat;
   q++; p++;
   }
   return result_sym;
  }
 } else {

 /* must be multi-D */
 /*several subscripts, the count (narg-1) can't be more than nd */
 if (narg > (nd+1) ) return execute_error(36);

 stride = 1;
 s_offset = 0;
 /* here we are manipulating symbol sets, so use pointers rather than
 the addresses used for real arrays, hence stride starts as 1
 we finally get the offsets for the symbol descriptors */
 for (k = 1; k < narg; k++) {			/*loop over arguments */
 iq = ps[k];	j = h->dims[k-1];
 switch (sym[iq].class)	{	/*switch on class of this subscript */

 case SCAL_PTR: iq = class8_to_1(iq);

 case SCALAR: 				/*scalar*/
  i = int_arg(iq);
  if ( i >= j || i < 0 ) return strarr_error(i, nsym);
  s_offset = s_offset + (i*stride);	stride = stride*j;
  break;

 case SUBSC_DESC: /*special subscript type (only allowed type here), we do not
	 return to the k loop, instead a new loop is started where the
	 old one left off, a return is executed before we go back */
  nd = nd - k + 1;		/* decrease nd by amount already done */
  /* this is similar to the array case but we are dealing with symbol pointers
  and compression is not supported (not sensible) */

  for (j=0;j<10;j++) { loop[j] = 0; rear[j] = -1; }		/*set ups */
  /* loop[0] = 0; */

  /* base[0] = (long) ptr; */
  base[0] = s_offset;	/* here as offsets from start, not absolute pointers */
  ncomp = nbtriv = nbdim = ind = 0;
  /* ind is dimension count in non-simple mode */
  /* nbtriv is the # of trivial dimensions */
  while (k < narg) {	/* loop over remaining dimensions */
  press[ind] = 0;
  inc[ind] = stride;	/* save stride in inc array for each dim. */
  j = h->dims[k-1];
  stride = stride*j;
  iq = ps[k];		/* current index, check if simple */
  if (sym[iq].class == 8) iq = class8_to_1(iq);
  if (sym[iq].class == 1) {  i = int_arg(iq);		/* simple case */
	  /* a simple one in the list */
	  if ( i >= j || i < 0 ) return strarr_error(i, nsym);
	  loop[ind] = lb[ind] = 1;	comf[ind] = 0;
	  offset[ind] = i*inc[ind]; base[0] += offset[ind];
	  if ( rear[ind] >= 0 ) return execute_error(37);
	  rear[ind] = ind;	nbtriv += 1;}
 else {			/*non-simple case, check if the allowed type */
 int	len;
 if (sym[iq].class != SUBSC_DESC) return execute_error(39);
 /* the class SUBSC_DESC symbol contains 4 integers, pull them out */
 pd = (int *) ((char *)sym[iq].spec.array.ptr + sizeof(struct ahead));
 i = *pd++;	/* first one is the starting index */ 
 if ( i >= j || i < 0 ) return strarr_error(i, nsym);
 len = *pd++;	/* second is the ending index */
 /* printf("start = %d, len = %d\n", i, len); */
 if ( len < 0 ) len = j - i; else  {
   if ( len >= j || len < i ) return strarr_error(i, nsym);
	 len = len - i + 1; }
 loop[ind] = lb[ind] = len;
 /* keep the compression info here in case we figure out a sensible thing
 to do with it */
 if (comf[ind] = *pd++) ncomp++;
 offset[ind] = i*inc[ind]; base[0] += offset[ind];
 /*now check the rearrange (fourth quantity) */
 if ( len == 1 || comf[ind] != 0 ) {	/*must be trivial, check it */
	  if ( rear[ind] >= 0 ) return execute_error(37);
	  rear[ind] = ind;	nbtriv += 1;}
 else {	len = *pd++;	/*re-use len variable for rearrange target */
	 if ( len < 0 ) len = nbdim;
	 len = len + nbtriv; if ( rear[len] >= 0 ) return execute_error(37);
	 rear[len] = ind;	nbdim += 1; }
 }
 ind += 1;
 k++; } 
 /* still in case SUBSC_DESC section, looping over subscripts finished */
 if ( ind != nd ) {
 printf("dimension count error, outer dimensions will be assumed 0\n");
 printf("nd, ind = %d, %d\n",nd, ind); nd = ind; }
 for (k=0;k<nd;k++) if (rear[k] < 0 ) return execute_error(37);
 /*trickle up base */
 for (k=1;k<=nd;k++) base[k] = base[0];
 /*rearrange everybody */
 (void) rearr(inc);	(void) rearr(loop);	(void) rearr(lb);	
 /*setup output array, get actual dimensions in scratch */
 j = 0;
 for (k=0;k<nd;k++) if (loop[k] > 1 && comf[k] == 0 ) scratch[j++]=loop[k];
 j = inc[0];		/*inner loop increment */
 if (ncomp) {
 /* compression (summation) not supported here */
	return execute_error(128); 
 } else {
  /* non-compression (always at present), if reduced to a single symbol,
  we just return the symbol rather than a symbol array */
  ns = 1; for (k=0;k<nbdim;k++) ns *= scratch[k];
  if (ns == 1) {
  	/* just get offset and break to single symbol stuff at end */
	s_offset = base[0];  break;
  } else {
  result_sym = symarr_scratch(nbdim, scratch);
  q = (struct sym_desc *)((char *)sym[result_sym].spec.array.ptr+sizeof(struct ahead));
  }
 
 while (1) { s_offset = base[1]; ns = lb[0];
 /* inner loop */
 while (ns) {
	stat = symbol_clone_via_ptr(q, p + s_offset, 1);
	if (stat != 1) return stat;
	q++;
 	s_offset += j; ns--; }
 i = 1;
 /*process outer counts */
	 while (1) { loop[i] -= 1;
	 if ( loop[i] > 0 ) break;
	 if (loop[i] < 0) return result_sym;
	 i++; }
 base[i] += inc[i];
 while (i > 1) { i--; base[i] = base [i+1]; loop[i] = lb[i]; }
 }					/*end of outer loop while */
 }
 }				/*end of switch */
 }
 /* we get here only for a series of scalar subscripts or a single symbol */
 /* check if offset in range */
 if (s_offset < 0 || s_offset >= nelem)  return execute_error(127);
 /* get the pointer location */
 /* printf("s_offset = %d\n", s_offset); */
 p =  p + s_offset;
 result_sym = find_next_temp();
 stat = symbol_clone_via_ptr(&sym[result_sym], p, 1); /* note mode = 1 here */
 if (stat != 1) return stat;
 return result_sym;
 }
 }
 /*------------------------------------------------------------------------- */
int ana_symclass(narg,ps)
 /*return the symbol class*/
 int narg,ps[];
 {
 int	nsym;
 nsym = ps[0];				/*the target symbol is the first */
 result_sym = scalar_scratch(2);
 sym[result_sym].spec.scalar.l = sym[nsym].class;
 return result_sym;
 }
 /*------------------------------------------------------------------------- */
int ana_isscalar(narg,ps)
 /*returns 1 if symbol is a scalar*/
 int narg,ps[];
 {
 int	nsym;
 nsym = ps[0];	 			/*the target symbol is the first */
 if ( sym[nsym].class==1 ) return 1; else return 4;
 }
 /*------------------------------------------------------------------------- */
int ana_isarray(narg,ps)
 /*returns 1 if symbol is a array*/
 int narg,ps[];
 {
 int	nsym;
 nsym = ps[0];	 			/*the target symbol is the first */
 if ( sym[nsym].class==4 ) return 1; else return 4;
 }
 /*------------------------------------------------------------------------- */
int ana_isstring(narg,ps)
 /*returns 1 if symbol is a string*/
 int narg,ps[];
 {
 int	nsym;
 nsym = ps[0];	 			/*the target symbol is the first */
 if ( sym[nsym].class==2 ) return 1; else return 4;
 }
 /*------------------------------------------------------------------------- */
int ana_isstrarr(narg,ps)
 /*returns 1 if symbol is a string array*/
 int narg,ps[];
 {
 int	nsym;
 nsym = ps[0];	 			/*the target symbol is the first */
 if ( sym[nsym].class==STRING_PTR ) return 1; else return 4;
 }
 /*------------------------------------------------------------------------- */
int ana_issymarr(narg,ps)
 /*returns 1 if symbol is a string array*/
 int narg,ps[];
 {
 int	nsym;
 nsym = ps[0];	 			/*the target symbol is the first */
 if ( sym[nsym].class==SYM_ARR ) return 1; else return 4;
 }
 /*------------------------------------------------------------------------- */
int ana_symdtype(narg,ps)
 /*return the symbol type*/
 int narg,ps[];
 {
 int	nsym;
 nsym = ps[0];				/*the target symbol is the first */
 result_sym = scalar_scratch(2);
 sym[result_sym].spec.scalar.l = sym[nsym].type;
 return result_sym;
 }
 /*------------------------------------------------------------------------- */
int ana_symaddress(narg,ps)
 /*return the symbol address*/
 int narg,ps[];
 {
 int	nsym;
 nsym = ps[0];				/*the target symbol is the first */
 result_sym = scalar_scratch(2);
 sym[result_sym].spec.scalar.l = (long) &sym[nsym].spec.scalar.l;
 return result_sym;
 }
 /*------------------------------------------------------------------------- */
int ana_num_elem(narg,ps)
 /*return the number of elements in the argument symbol */
 int narg,ps[];
 {
 int	nsym, nd, n, j;
 struct	ahead	*h;
 nsym = ps[0];				/*the target symbol is the first */
 result_sym = scalar_scratch(2);
 switch (sym[nsym].class)	{			/*switch on class */
 case STRING_PTR:	/* works just like array */
 case ARRAY:
 case SYM_ARR:
	 h = (struct ahead *) sym[nsym].spec.array.ptr;
	 nd = h->ndim;
	 n = 1; for (j=0;j<nd;j++) n *= h->dims[j]; /*# of elements for nsym */
	 break;
 case STRING:
	 n = sym[nsym].spec.array.bstore - 1;	/*don't include the null */
	 break;
 case SCAL_PTR:
 case SCALAR: n = 1;				/*scalar, just return 1 */
	 break;
 default: return execute_error(32);
 }
 sym[result_sym].spec.scalar.l = n;
 return result_sym;
 }
 /*----------------------------------------------------------------------- */
int ana_num_dimen(narg,ps)
 /*return the number of dimensions in the argument symbol */
 int narg,ps[];
 {
 int	nsym, nd;
 struct	ahead	*h;
 nsym = ps[0];				/*the target symbol is the first */
 result_sym = scalar_scratch(2);
 switch (sym[nsym].class)	{			/*switch on class */
 case STRING_PTR:	/* works just like array */
 case SYM_ARR:		/* works just like array */
 case ARRAY:
	 h = (struct ahead *) sym[nsym].spec.array.ptr;
	 nd = h->ndim;
	 break;
 case STRING:
	 nd = 1;					/*return 1 for strings */
	 break;
 case SCAL_PTR:
 case SCALAR: nd = 0;				/*scalar, just return 0 */
	 break;
 default: return execute_error(32);
 }
 sym[result_sym].spec.scalar.l = nd;
 return result_sym;
 }
 /*------------------------------------------------------------------------- */
int ana_dimen(narg,ps)
 /* return a particular dimension in the argument symbol */
 /* returns all dimensions if only a single argument is specified. */
 /* -- LS 26jan00 */
 int narg,ps[];
 {
 int	nsym, nd, n, i;
 struct	ahead	*h, *hout;
 nsym = ps[0];				/*the target symbol is the first */
 if (narg > 1) {
   i = int_arg(ps[1]);
   result_sym = scalar_scratch(2);
   switch (sym[nsym].class)	{			/*switch on class */
     case STRING_PTR:	/* works just like array */
     case ARRAY:
     case SYM_ARR:
       h = (struct ahead *) sym[nsym].spec.array.ptr;
       nd = h->ndim;
       if ( i >= 0 && i < nd ) n = h->dims[i]; else n = 0;
       break;
     case STRING:		/*string */
       if ( i == 0 )
	 n = sym[nsym].spec.array.bstore - 1;	/*don't include the null */
       else n = 0;
       break;
     case SCAL_PTR:
     case SCALAR: n = 0;	/*scalar, just return 0 */
       break;
     default: return execute_error(32);
   }
 } else {
   switch (sym[nsym].class)	{			/*switch on class */
     case STRING_PTR:	/* works just like array */
     case ARRAY:
     case SYM_ARR:
       h = (struct ahead *) sym[nsym].spec.array.ptr;
       nd = h->ndim;
       result_sym = array_scratch(2, 1, &nd);
       hout = (struct ahead *) sym[result_sym].spec.array.ptr;
       memcpy(hout + 1, h->dims, nd*sizeof(int));
       return result_sym;
       break;
     case STRING:		/*string */
       if ( i == 0 )
	 n = sym[nsym].spec.array.bstore - 1;	/*don't include the null */
       else n = 0;
       break;
     case SCAL_PTR:
     case SCALAR: n = 0;	/*scalar, just return 0 */
       break;
     default: return execute_error(32);
   }
 }
 sym[result_sym].spec.scalar.l = n;
 return result_sym;
 }
 /*------------------------------------------------------------------------- */
int ana_dimen_symarr(narg,ps)
 /* specialized routine to return a dimension of a symarr element without
 extracting it. The call nx = dimen_symarr(a, k, i) is logically the same
 as  nx = dimen(a(k), i) but the latter causes a copy of the data in a(i)
 to be made, this is undesirable if a(i) is a large image or worse.
 3 args are required */
 int narg,ps[];
 {
 int	nsym, nd, n, i, k, j;
 struct	ahead	*h, *hout;
 struct sym_desc *pdesc;
 nsym = ps[0];				/*the target symbol is the first */
 if (int_arg_stat(ps[1], &k) != 1) return -1;
 if (int_arg_stat(ps[2], &i) != 1) return -1;
 result_sym = scalar_scratch(2);
 switch (sym[nsym].class)	{			/*switch on class */
   case SYM_ARR:
     h = (struct ahead *) sym[nsym].spec.array.ptr;
     /* check the # of symbols against k */
     nd = h->ndim;
     n = 1;
     for(j=0;j<nd;j++) n *= h->dims[j];
     if (k >= n || k < 0) {
     	printf("out of range symarr index = %d\n", k);
     	return -1; }
     pdesc = (struct sym_desc *) ((char *)h + sizeof(struct ahead));
     pdesc += k;
     /* we must have an array here or taking the dimension will cause
     a problem */
     switch (pdesc->class)	{			/*switch on class */
       case ARRAY:
       case SYM_ARR:
	 h = (struct ahead *) pdesc->spec.array.ptr;
	 nd = h->ndim;
	 if ( i >= 0 && i < nd ) n = h->dims[i]; else n = 0;
	 break;
       default:
         n = 0;
     }
     break;
   default: return execute_error(32);
   }
 sym[result_sym].spec.scalar.l = n;
 return result_sym;
 }
 /*------------------------------------------------------------------------- */
int ana_redim_f(narg,ps)
 /* redimension an array, a function version, returns input symbol */
 int narg,ps[];
 {
 int	nsym;
 nsym = ps[0];				/*the target symbol is the first */
 if ( ana_redim(narg,ps) == 1 ) return nsym; else return -1;
 }
 /*------------------------------------------------------------------------- */
int ana_redim(narg,ps)
 /* redimension an array */
 int narg,ps[];
 {
 extern	int	redim_warn_flag;
 int	nsym, nd, n, j, mq, nn, pd[8], class;
 struct	ahead	*h;
 nsym = ps[0];				/*the target symbol is the first */
 /* 3/21/00 upgrade to handle strarr and symarr properly */
 class = sym[nsym].class;
 /* check if an array of some sort */
 switch (class) {
  case ARRAY:
  case STRING_PTR:
  case SYM_ARR:
   break;
  default: return execute_error(67);
 }
 n = nn = 1;
 /* get old size */
 h = (struct ahead *) sym[nsym].spec.array.ptr;
 nd = h->ndim;
 for (j=0;j<nd;j++) { nn *= h->dims[j];  } 
 /* nn is now the # of original elements */

 if ( (nd = narg - 1) <= 0 ) return execute_error(62);
 if ( nd > 8 ) return execute_error(63);
 if ( get_dims(&nd,&ps[1],pd) != 1) return -1;
 for (j=0;j<nd;j++) { n *= pd[j]; }			/*get the total */
 /* n is now the # of new elements */

 switch (class) {
  case ARRAY:
   j = ana_type_size[ sym[nsym].type ];
   n = n*j;  nn = nn*j;
   /* compare with storage for array rather than current dimensions in case it
   was already changed but data is still there */
   mq = sym[nsym].spec.array.bstore - sizeof( struct ahead );
   if ( n > mq ) return execute_error(64);
   break;
  case STRING_PTR:
  case SYM_ARR:
   /* here we must always be the same or smaller, also we have to be
   irreversible and delete any strings or symbols that are no longer
   in the set */
   if (n > nn) return execute_error(64);
   if (n < nn) {
    /* OK, get rid of the unwanted ones at the end */
    int nelem;
    nelem = nn - n;
    if (class == STRING_PTR) {
    char	**p;
    p = (char **) ((char *)h + sizeof(struct ahead));
    p = p + n;
    while (nelem--) { if (*p) free(*p); p++; }
    } else

    if (class == SYM_ARR) {
    struct sym_desc *p;
    int	stat;
    p = (struct sym_desc *) ((char *)h + sizeof(struct ahead));
    p = p + n;
    while (nelem--) { 
     /* each symbol in the array must be deleted */
     stat = delete_symbol_via_ptr(p);
     if (stat != 1) return stat;
     p++;
     }
     }
    }
   break;
  default: return execute_error(67);
  }

 /* proceed, also give warning if total now smaller then old dimensions */
 h->ndim = nd;
 for (j=0;j<nd;j++) { h->dims[j] = pd[j]; }
 if (n < nn && redim_warn_flag)
	 printf("WARNING: result from REDIM is smaller than original\n");
 return 1;
 }
 /*----------------------------------------------------------------------- */
int ana_concat(narg,ps)
 /* 3/26/97 strings are now concatenated into string arrays */
 /* ana's concatenation function, insisted on by idl users, but works a bit
	 differently
 there are 2 modes
 mode 1 - has a single argument, we return a variable which has one extra
	  dimension which has a length of one; i.e., a scalar becomes a
	  vector of length 1, a vector of length n becomes an array of nx1
 mode 2 - 2 or more arguments,result combines the last dimension
	 all dimensions but the last must be of the same size */
 /* C version 11/2/91, r shine */
 int narg,ps[];
 {
 union	types_ptr q1,q2;
 char	**p, **q, *p2, *p1;
 int	nd, j, i, dim[8];
 int	iq, nsym, mq=0, toptype=0, topnd=0, sflag=0, n, nq;
 struct	ahead	*h;
 if ( narg <=0 ) return execute_error(55);
 if ( narg == 1)	{				/* one argument case */
   if ( (iq = ps[0]) <= 0 ) return iq;			/* error check */
   switch (sym[iq].class)	{
   case 8:						/*scalar ptr */
    iq = class8_to_1(iq);			/*just eval it and treat as scalar */
   case 1:						/*scalar case */
    /* create an array of length 1 */
    dim[0]=1;	nsym = array_scratch(sym[iq].type,1,dim);
     bcopy((char *) &sym[iq].spec.scalar.l,
	   (char *) ((char *) sym[nsym].spec.array.ptr + sizeof(struct ahead)),
	   ana_type_size[sym[iq].type]);
    return nsym;
   case STRING:						/*string */
    /* make a string array with just one element */
    dim[0]=1;	nsym = strarr_scratch(1, dim);
    p = (char **) ((char *) sym[nsym].spec.array.ptr + sizeof(struct ahead));
    p2 = *p;	/* just the first and only element */
    p1 = (char *) sym[iq].spec.array.ptr;
    *p = strsave(p1);
    return nsym;
   case STRING_PTR:					/*strarr */
    { char **p, **q;
    /* make a copy with an additional dimension (of size 1) added */
    h = (struct ahead *) sym[iq].spec.array.ptr;
    q = (char **) ((char *) h + sizeof(struct ahead));
    nd = h->ndim;
    n = 1; for (j=0;j<nd;j++) { dim[j] = h->dims[j]; n *= h->dims[j]; }
    dim[nd] = 1;	nd += 1;
    nsym = strarr_scratch(nd, dim);
    h = (struct ahead *) sym[nsym].spec.array.ptr;
    /* now the destination */
    p = (char **) ((char *) sym[iq].spec.array.ptr + sizeof(struct ahead));
    mq = sym[iq].spec.array.bstore - sizeof(struct ahead);
    bcopy( (char *)sym[iq].spec.array.ptr + sizeof(struct ahead),
  	  (char *)sym[nsym].spec.array.ptr + sizeof(struct ahead), mq);
    /* and copy the strings themselves */
    while (n--) { if (*q) *p = strsave(*q);  p++;  q++; }
    return nsym;
   case ARRAY:						/*array */
    /* make a copy with an additional dimension (of size 1) added */
    h = (struct ahead *) sym[iq].spec.array.ptr;
    q1.l = (int *) ((char *)h + sizeof(struct ahead));
    nd = h->ndim;
    for (j=0;j<nd;j++) dim[j] = h->dims[j];
    dim[nd] = 1;	nd += 1;	nsym = array_scratch(sym[iq].type,nd,dim);
    h = (struct ahead *) sym[nsym].spec.array.ptr;
    q2.l = (int *) ((char *)h + sizeof(struct ahead));
    /* 8/14/2000 - the following had the bstore for iq rather than nsym, this
    is wrong if iq was truncated via a redim. Caused a hard to find bug! */
    mq = sym[nsym].spec.array.bstore - sizeof(struct ahead);
    bcopy( (char *) q1.l, (char *) q2.l, mq);
    return nsym;
    }
   default: return execute_error(57);
   }						/*end of switch */

 } else {					/* 2 or more arguments */

 /* first a quick look at all the symbols */
 for (i=0;i<narg;i++) {
 if ( (iq = ps[i]) <= 0 ) return iq;			/* error check */
 switch (sym[iq].class)	{
 case 8:						/*scalar ptr */
 case 1:						/*scalar case */
  mq += 1;	break;
 case STRING:					/*string */
  /* 3/26/97 strings are now combined into a string array */
  mq += 1;
  sflag += 1;					/* a string flag */
  break;
 case STRING_PTR:				/*strarr */
  sflag += 1;					/* a string flag */
  /* rest the same as an array for this pass */
 case ARRAY:					/*array */
  h = (struct ahead *) sym[iq].spec.array.ptr;
  nd = h->ndim;
  n = 1; for (j=0;j<nd;j++) n *= h->dims[j]; /*# of elements for array */
  if ( nd > topnd ) topnd = nd;
  mq += n;
  break;
 default: return execute_error(57);
 }
 if (sym[iq].type > toptype ) toptype = sym[iq].type;
 }

 /* now run through again */

 if (sflag )	{			/*string case, all must be strings */
  if ( sflag < narg ) return execute_error(58);	/*not all were */
  /* similar to arrays (below) with strings as scalars and strarr's as arrays*/
  /* use first one to set up inner dimensions */
  iq = ps[0];	n = MAX(0, topnd-1);
  if (sym[iq].class == 2) {	/* first a string */
   /* must be all string or a string/strarr combination */
   if (topnd > 1) return execute_error(59);
   dim[0] = mq; } else {			/* first a strarr */
   if (sym[iq].class != STRING_PTR) return execute_error(57);
   h = (struct ahead *) sym[iq].spec.array.ptr;
   nd = h->ndim;
   if ( nd != topnd && nd != n) return execute_error(59);
   nq = 1; for (j=0;j<n;j++) { dim[j] = h->dims[j];  nq *= h->dims[j]; }
   if ( mq%nq ) return execute_error(59);
   dim[n] = mq/nq; }
  /* each argument must have a dimension=topnd or topnd-1 and inner n must
	  match, check before allocating memory */
  for (i=0;i<narg;i++) {
   iq = ps[i];
   if (sym[iq].class == STRING_PTR) {
    h = (struct ahead *) sym[iq].spec.array.ptr;	nd = h->ndim; } else {
    nd = 0; }
   if ( nd != topnd && nd != topnd-1) return execute_error(59);
   if (n) for (j=0;j<n;j++) if (dim[j] != h->dims[j]) return execute_error(59);
  }
  /* get the resultant array */
  nsym = strarr_scratch(MAX(topnd,1), dim);
  h = (struct ahead *) sym[nsym].spec.array.ptr;
  p = (char **) ((char *) h + sizeof(struct ahead));
  /* now just load everybody in */
  for (i=0;i<narg;i++) {
   iq = ps[i];
   switch (sym[iq].class)	{
   case STRING:
    p1 = (char *) sym[iq].spec.array.ptr;
    *p++ = strsave(p1);
    break;
   case STRING_PTR:
    /* grab each one and copy */
    h = (struct ahead *) sym[iq].spec.array.ptr;
    nd = h->ndim;
    n = 1; for (j=0;j<nd;j++) n *= h->dims[j];
    q = (char **) ((char *) h + sizeof(struct ahead));
    while (n--) { if (*q) *p = strsave(*q);  p++;  q++; }
    break;
   default: return execute_error(58); /*just in case */
   }
  }
  return	nsym;
  }


 if ( topnd == 0 ) {				/* all scalars case */
  dim[0] = mq;	nd = 1;	nsym = array_scratch(toptype,nd,dim);
  n = ana_type_size[toptype];
  h = (struct ahead *) sym[nsym].spec.array.ptr;
  q2.b = (byte*) ((char *)h + sizeof(struct ahead));
  for (i=0;i<narg;i++) {
   iq = ps[i];
   if ( sym[iq].class == 8 ) iq = class8_to_1(iq);
   if ( sym[iq].class != 1 ) return execute_error(57); /*just in case */
   bcopy( (char *) &sym[iq].spec.scalar.l, (char *) q2.b, n);
   if (sym[iq].type != toptype ) ana_convert( q2.b, sym[iq].type, toptype);
   q2.b += n;		/* add the right byte count for each element */
 }
 return	nsym;
 }
 /* arrays, the only possibility left */
 /* use first one to set up inner dimensions */
 iq = ps[0];	n = topnd-1;
 if (sym[iq].class == 8 || sym[iq].class == 1) {	/* first a scalar */
  /* must be a scalar/vector combination */
  if (topnd != 1) return execute_error(59);
  dim[0] = mq; } else {				/* first an array */
  if (sym[iq].class != 4) return execute_error(57);
  h = (struct ahead *) sym[iq].spec.array.ptr;
  nd = h->ndim;
  if ( nd != topnd && nd != n) return execute_error(59);
  nq = 1; for (j=0;j<n;j++) { dim[j] = h->dims[j];  nq *= h->dims[j]; }
  if ( mq%nq ) return execute_error(59);
  dim[n] = mq/nq; }
 /* each argument must have a dimension=topnd or topnd-1 and inner n must
	 match, check before allocating memory */
 for (i=0;i<narg;i++) {
  iq = ps[i];
  if (sym[iq].class == 4) {
   h = (struct ahead *) sym[iq].spec.array.ptr;	nd = h->ndim; } else {
   nd = 0; }
  if ( nd != topnd && nd != topnd-1) return execute_error(59);
  if (n) for (j=0;j<n;j++) if (dim[j] != h->dims[j]) return execute_error(59);
 }
 /* get the resultant array */
 nsym = array_scratch(toptype, topnd, dim);
 h = (struct ahead *) sym[nsym].spec.array.ptr;
 q2.b = (byte*) ((char *)h + sizeof(struct ahead));
 /* now just load everybody in, converting type as necessary */
 for (i=0;i<narg;i++) {
  iq = ps[i];
  switch (sym[iq].class)	{
   case 8:						/*scalar ptr */
    iq = class8_to_1(iq);
   case 1:						/*scalar case */
    n = 1;	q1.l = &sym[iq].spec.scalar.l;	break;
   case 4:						/*array */
    h = (struct ahead *) sym[iq].spec.array.ptr;
    q1.l = (int *) ((char *)h + sizeof(struct ahead));
    nd = h->ndim;
    n = 1; for (j=0;j<nd;j++) n *= h->dims[j];	break;
   }
  switch (toptype) {
   case 0: while (n) { *q2.b++ =  (*q1.b++);n--;} break;
   case 1: switch (sym[iq].type) {
	 case 1: while (n) { *q2.w++ =  (*q1.w++);n--;} break;
	 case 0: while (n) { *q2.w++ = (short) (*q1.b++);n--;} break;
	 }	break;
   case 2: switch (sym[iq].type) {
	 case 2: while (n) { *q2.l++ =  (*q1.l++);n--;} break;
	 case 1: while (n) { *q2.l++ = (int) (*q1.w++);n--;} break;
	 case 0: while (n) { *q2.l++ = (int) (*q1.b++);n--;} break;
	 }	break;
   case 3: switch (sym[iq].type) {
	 case 3: while (n) { *q2.f++ =  (*q1.f++);n--;} break;
	 case 2: while (n) { *q2.f++ = (float) (*q1.l++);n--;} break;
	 case 1: while (n) { *q2.f++ = (float) (*q1.w++);n--;} break;
	 case 0: while (n) { *q2.f++ = (float) (*q1.b++);n--;} break;
	 }	break;
   case 4: switch (sym[iq].type) {
	 case 4: while (n) { *q2.d++ =  (*q1.d++);n--;} break;
	 case 3: while (n) { *q2.d++ = (double) (*q1.f++);n--;} break;
	 case 2: while (n) { *q2.d++ = (double) (*q1.l++);n--;} break;
	 case 1: while (n) { *q2.d++ = (double) (*q1.w++);n--;} break;
	 case 0: while (n) { *q2.d++ = (double) (*q1.b++);n--;} break;
	 }	break;
  }
 }
 return nsym;
 }
 }						/*end of ana_concat */
 /*------------------------------------------------------------------------- */
int ana_symarr_merge(narg,ps)
 int narg,ps[];
 /* a way to combine two symbol arrays without copying their contents,
 the call is x = symarr_merge(a,b) where a and b must be symbol arrays,
 x becomes a concatenation of the two, copying their descriptor arrays
 with those of b attached to the end of the a list. The pointers to the
 target symbols are intact. Then a and b are destroyed! This prevents
 any questions about who can redefine the targets. Not very elegant. The
 complement operation is symarr_split which is coded as a subroutine */
 {
 int	iq, jq, *p1, *p2, n1, n2, dim[8], nd, result_sym, n;
 struct	ahead	*h;
 struct sym_desc *pd1, *pd2, *pdout;
 iq = ps[0];
 jq = ps[1];
 if (sym[iq].class != SYM_ARR || sym[jq].class != SYM_ARR) {
 	printf("arguments must be symbol arrays\n");
 	return -1; }
 /* a limitation here is that the result will be a 1-D symbol array
 regardless of what a and b are, this could be extended to work like
 regular arrays with some extra effort. All inner dimensions would have
 to agree though. */
 n1 = get_p_c(&iq, &p1);
 n2 = get_p_c(&jq, &p2);
 pd1 = (struct sym_desc *) p1;
 pd2 = (struct sym_desc *) p2;
 /* get the output symbol */
 nd = 1;
 n = n1 + n2;
 dim[0] = n;
 result_sym = symarr_scratch(nd, dim);
 pdout = (struct sym_desc *)((char *)sym[result_sym].spec.array.ptr+sizeof(struct ahead));
 /* devour the descriptors */
 while (n1--) { *pdout++ = *pd1;  pd1->class = 0; pd1++; }
 while (n2--) { *pdout++ = *pd2;  pd2->class = 0; pd2++; }
 /* and toss the bones (delete the input symbols) */
 if ( delete_symbol(iq) != 1 ) return execute_error(17);
 if ( delete_symbol(jq) != 1 ) return execute_error(17);
 return result_sym;
 }
 /*------------------------------------------------------------------------- */
int ana_symarr_split(narg,ps)
 int narg,ps[];
 /* a subroutine */
 /* a way to split a symbol array into 2 new ones without copying the contents,
 the call is symarr_split, x, a, b, n where x must be a symbol array of
 length n+1 or greater. x gets destroyed and a gets the first n elements
 while b gets whatever is left.
 */
 {
 int	sq, iq, jq, *p0, n0, n1, n2, dim[8], nd, n;
 int	flag1=0, flag2=0;
 struct	ahead	*h;
 struct sym_desc *pd1, *pd2, *pdin;
 sq = ps[0];
 iq = ps[1];
 jq = ps[2];
 printf("sq,iq,jq = %d %d %d\n", sq,iq,jq);
 if (iq == jq) {
 	printf("the 2 destination symbols cannot be the same\n");
 	return -1; }
 if (sym[sq].class != SYM_ARR) {
 	printf("argument must be a symbol array\n");
 	return -1; }
 /* a limitation here is that the results will be 1-D symbol arrays
 regardless of what x was, this could be extended to work like
 regular arrays with some extra effort. */
 n0 = get_p_c(&sq, &p0);
 pdin = (struct sym_desc *) p0;
 if (int_arg_stat(ps[3], &n) != 1) return -1;
 if (n >= n0 || n < 0 ) {
 	printf("position %d too large, must be <= %d\n", n, n0);
 	return -1; }
 /* create (redefine) the sym arrays a and b, alias iq and jq here */
 /* need to check if one is the same as "input" symbol, in this case we
 have to create a temp. We don't do this otherwise in order to free any
 memory no longer needed before allocating more. */
 nd = 1;
 n1 = dim[0] = n;
 if (iq != sq) {
  if ( redef_symarr(iq, nd, dim) != 1 ) return -1;
  } else {
  flag1 = 1;
  iq = symarr_scratch(nd, dim);
  }
 pd1 = (struct sym_desc *)((char *)sym[iq].spec.array.ptr+sizeof(struct ahead));
 nd = 1;
 n2 = dim[0] = n0 - n;
 if (jq != sq) {
 if ( redef_symarr(jq, nd, dim) != 1 ) return -1;
  } else {
  flag2 = 1;
  jq = symarr_scratch(nd, dim);
  }
 pd2 = (struct sym_desc *)((char *)sym[jq].spec.array.ptr+sizeof(struct ahead));
 while (n1--) { *pd1++ = *pdin;  pdin->class = 0; pdin++; }
 while (n2--) { *pd2++ = *pdin;  pdin->class = 0; pdin++; }

 /* if one of the results was the same symbol as the input, do the equate
 now, else we just delete the input */
 //printf("flag1, flag2 = %d, %d\n", flag1, flag2);
 if (flag1) { if (ana_replace(sq, iq) != 1) return -1;
  } else if (flag2) { if (ana_replace(sq, jq) != 1) return -1;
   } else if ( delete_symbol(sq) != 1 ) return execute_error(17);
 return 1;
 }
 /*------------------------------------------------------------------------- */
int ana_symarr_mask(narg,ps)
 int narg,ps[];
 /* a subroutine */
 /* eliminates entries in a symbol array without copying,
 the call is symarr_mask, x, ind where x must be a symbol array and ind
 is a mask array of the indices to keep in the modified version of x, the
 entries equal to zero in ind are eliminated and their memory allocations
 recovered. The symbol array pointer array is also recreated and the
 original one deleted. Note that ind must match x in # of elements.
 */
 {
 int	sq, iq, jq, *p0, *p1, *p, n0, n1, n2, dim[8], nd, n;
 int	flag1=0, flag2=0;
 struct	ahead	*h;
 struct sym_desc *pd1, *pdin;
 sq = ps[0];
 ind = ps[1];
 if (sym[sq].class != SYM_ARR) {
 	printf("argument must be a symbol array\n");
 	return -1; }
 /* a limitation here is that the result will be a 1-D symbol array
 regardless of what x was */
 n0 = get_p_c(&sq, &p0);
 pdin = (struct sym_desc *) p0;
 /* check ind and its size */
 if (sym[ind].class != ARRAY) return execute_error(66);
 ind = ana_long(1, &ps[1] );
 n1 = get_p_c(&ind, &p1);
 /* the two must match */
 if (n1 != n0) {
   printf("size mismatch between symbol array and mask: %d vs %d\n",n0,n1);
   return -1; }
 
 /* determine how many we need, if all are selected we just return (nothing
 to do */
 n = 0;
 p = p1;
 while (n1--) if (*p++) n++;
 if (n == n0) return 1;
 
 /* create the modified sym array x, alias iq here */
 nd = 1;
 n1 = dim[0] = n;
 iq = symarr_scratch(nd, dim);

 pd1 = (struct sym_desc *)((char *)sym[iq].spec.array.ptr+sizeof(struct ahead));

 /* the ones we want to save are copied to the new symbol array and the
 class is zeroed in the original so the symbol is "hidden" when we delete
 the original symbol array and all its non-hidden symbols */
 p = p1;
 while (n0--) {
 	if (*p++) { *pd1++ = *pdin; pdin->class = 0; }
 	pdin++; }

 if (ana_replace(sq, iq) != 1) return -1;
 return 1;
 }
 /*------------------------------------------------------------------------- */