/* file subsc.c  ana subscript routines */
#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 */
#ifdef __alpha
 static long	base[10];	/* long is 64 bits on alpha */
#ifdef __sgi
#if _MIPS_SZPTR == 64
 static	long	base[10];
#endif
#endif
#else
 static int	base[10];
#endif
 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, nx, ny, nd1, nd2, type1, type2, ix, iy, nx2, ny2, nc, del, del2;
 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 = 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 */
 iq = (iy * nx + ix);
 del2 = MAX( 0, nx2 - nx + ix);		/* needed if overrun */
 nx2 = MIN( nx - ix, nx2);		/* transfer count per row */
 del = nx - nx2;				/*leftovers per row */
 ny2 = MIN( ny2, ny - iy);	/* was a bug here 11/25/92 */
 if (ny2 < 1 || nx2 < 1) { printf("warning - insert out of bounds\n");
		 return 1; }
 /* use only 4 switches for 4 sizes */
 switch (type1) {
 case 0:
 p1.b += iq;
 while (ny2--) { nc = nx2; 
   while (nc--) *p1.b++ = *p2.b++;  p1.b += del; p2.b += del2; } break;
 case 1:
 p1.w += iq;
 while (ny2--) { nc = nx2; 
   while (nc--) *p1.w++ = *p2.w++;  p1.w += del; p2.w += del2; } break;
 case 2:
 case 3:
 p1.l += iq;
 while (ny2--) { nc = nx2; 
   while (nc--) *p1.l++ = *p2.l++;  p1.l += del; p2.l += del2; } break;
 case 4:
 p1.d += iq;
 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];

  /* 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); break;
 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 */
	/* simple for initial configuration, just allow a single subscript
	and return a copy of the string */
	return strptr_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 9:					/*symbol ptr. 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 */
  dim[0]=ns;	result_sym = strarr_scratch(1, dim);
  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 = n - i; else  {
	 if ( ns >= n || 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;	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] = 0;		/* 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));
  }
 
 
 /* no compression case */
 /* switch on source type */
 /* outer loops */
 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 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_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:
	 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 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 */
 int narg,ps[];
 {
 int	nsym, nd, n, i;
 struct	ahead	*h;
 nsym = ps[0];				/*the target symbol is the first */
 if (narg > 1) i = int_arg(ps[1]); else i = 0; /*default is first dimension */
 result_sym = scalar_scratch(2);
 switch (sym[nsym].class)	{			/*switch on class */
 case STRING_PTR:	/* works just like array */
 case ARRAY:
	 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); break;
 }
 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];
 struct	ahead	*h;
 nsym = ps[0];				/*the target symbol is the first */
 if (sym[nsym].class != ARRAY &&  sym[nsym].class != STRING_PTR)
 	return execute_error(67);
 /* get old size */
 nn = n = ana_type_size[ sym[nsym].type ];
 h = (struct ahead *) sym[nsym].spec.array.ptr;
 nd = h->ndim;
 for (j=0;j<nd;j++) { nn *= h->dims[j];  } 
 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 */
 /* 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);
 /* 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[10];
 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));
  mq = sym[iq].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 */
 /*------------------------------------------------------------------------- */