/*a collection of internal ana subroutines and functions */
#include <stdio.h>
#include "ana_structures.h"
#include "defs.h"
#include <time.h>
#include <stdlib.h>
#include <sys/time.h>
#include <string.h>
#include <math.h>
/* a better conditional than APPLE might be available */
#if __APPLE__
#include <Accelerate/Accelerate.h>
#endif

 extern	struct sym_desc sym[];
 extern	int	temp_base, n_symbols;
 extern	char *find_sym_name(int iq);
 extern	char	*strsave();
 extern	int	vfix_top;
 extern	float	float_arg();
 extern	int	ana_type_size[];
 union	types_ptr { byte *b; short *w; int *l; float *f; double *d;};		
 static	int	result_sym;
 static	float	scount;
 double	(*mc[])() = { sin, cos, tan, asin, acos, atan, sinh, cosh, tanh,
	 sqrt, cbrt, exp, expm1, log, log10, log1p,
	 erf, erfc, atan2, j0, j1, y0, y1, pow, yn, jn, hypot};
 extern	therm2temp(), patan();
 int	(*mcint[])() = { therm2temp, patan };
 int	num_of_nans;
 static char	*month_names[12] = {"Jan","Feb","Mar","Apr","May","Jun",
 		"Jul","Aug","Sep","Oct","Nov","Dec" };
 /*------------------------------------------------------------------------- */
int get_p_c_r(int *nsym, int **p, int *result_sym, int **presult)
 /* returns the number of elements for a symbol and a pointer to the
 first one, returns 0 if something is wrong; will change the symbol
 number if a class 8
 also returns a result symbol which is the same type and size as nsym,
 if nsym is a temp, it will be exactly the same (reuse) - hence be careful */
 {
 int	n, nd, i;
 struct	ahead	*h;
 switch (sym[*nsym].class)	{
 /* following note copied from get_p_c, implemented here also but
 doesn't change much for this utility since the input symbol is
 generally not considered changable in routines that call here */
 /* 5/11/96 the class 8 used to call class8_to_1 which made a copy
 of the class8. This is fine for input only. To allow the value to be
 changed, we now pass the actual class pointer. This also means that
 the symbol is not changed but keep the call by reference in case we
 need it to handle new situations. If that doesn't occur, consider
 changing to a value call. */
 case SCAL_PTR:  					/*scalar ptr case */
 	*p = (int *)sym[*nsym].spec.array.ptr;
	*result_sym = scalar_scratch(sym[*nsym].type);
 	*presult = &sym[*result_sym].spec.scalar.l;
	return 1;
 case SCALAR:					/*scalar case */
	if ( *nsym >= temp_base ) { *result_sym = *nsym; }
	else *result_sym = scalar_scratch(sym[*nsym].type);
 	*p = &sym[*nsym].spec.scalar.l;
 	*presult = &sym[*result_sym].spec.scalar.l;
	return 1;
 case ARRAY:
	h = (struct ahead *) sym[*nsym].spec.array.ptr;
	*p = (int *) ((char *)h + sizeof(struct ahead));
	nd = h->ndim;
	n = 1; for (i=0;i<nd;i++) n *= h->dims[i]; /*# of elements for nsym */
	if ( *nsym >= temp_base ) { *result_sym = *nsym; }
	else *result_sym = array_clone(*nsym, sym[*nsym].type);
	h = (struct ahead *) sym[*result_sym].spec.array.ptr;
	*presult = (int *) ((char *)h + sizeof(struct ahead));
	return n;
 case STRING:
 	return execute_error(47);
 default:
 	return execute_error(48);
 }
 }
 /*------------------------------------------------------------------------- */
int get_p_c(int *nsym, int **p)
 /* returns the number of elements for a symbol and a pointer to the
 first one, returns 0 if something is wrong */
 {
 int	n, nd, i;
 struct	ahead	*h;
 switch (sym[*nsym].class)	{
 /* 5/11/96 the class 8 used to call class8_to_1 which made a copy
 of the class8. This is fine for input only. To allow the value to be
 changed, we now pass the actual class pointer. This also means that
 the symbol is not changed but keep the call by reference in case we
 need it to handle new situations. If that doesn't occur, consider
 changing to a value call. */
 case SCAL_PTR:  					/*scalar ptr case */
 	*p = (int *)sym[*nsym].spec.array.ptr;
	return 1;
 case SCALAR:					/*scalar case */
 	*p = &sym[*nsym].spec.scalar.l;
	return 1;
 case ARRAY:
 case SYM_ARR:		/* works just like array */
	h = (struct ahead *) sym[*nsym].spec.array.ptr;
	*p = (int *) ((char *)h + sizeof(struct ahead));
	nd = h->ndim;
	n = 1; for (i=0;i<nd;i++) n *= h->dims[i]; /*# of elements for nsym */
	return n;
 case STRING:
 	return execute_error(47);
 default:
 	return execute_error(48);
 }
 }
 /*------------------------------------------------------------------------- */
int ana_defined(narg,ps)
 /* returns 1 if defined, 0 if not */
 int narg, ps[];
 {
 int	iq;
 iq = ps[0];
 if ( sym[iq].class == 255 || sym[iq].class == 0) return 4; else return 1;
 }
 /*------------------------------------------------------------------------- */
int ana_notdefined(narg,ps)
 /* returns 0 if defined, 1 if not */
 int narg, ps[];
 {
 int	iq;
 iq = ps[0];
 if ( sym[iq].class == 255 || sym[iq].class == 0) return 1; else return 4;
 }
 /*------------------------------------------------------------------------- */
int ana_zeroifnotdefined(narg,ps)
 /* if symbol in list is not defined, set to a scalar zero*/
 int narg, ps[];
 {
 int iq, i, z=0;
 char *p;
 for (i=0;i<narg;i++) {
   iq = ps[i];
   if ( sym[iq].class == 255 || sym[iq].class == 0) {
     /* not defined, so do it */
     if (redef_scalar( iq, 2, &z) != 1) return execute_error(13);
   }
 }
 return 1;
 }
 /*------------------------------------------------------------------------- */
int ana_zapsym(narg,ps)
 /* defines all symbols in list to a scalar zero, used to recover memory */
 int narg, ps[];
 {
 int iq, i, z=0;
 char *p;
 for (i=0;i<narg;i++) {
   iq = ps[i];
   if (redef_scalar( iq, 2, &z) != 1) return execute_error(13);
 }
 return 1;
 }
 /*------------------------------------------------------------------------- */
int ana_quit(narg,ps)
 /*exit routine, calls are exit,status or quit,status */
 int narg, ps[];
 {
 int	iq;
 if (narg) iq = int_arg(ps[0]); else iq = 1;
 (void) ana_exit(iq);
 return 0;
 }
 /*------------------------------------------------------------------------- */
int ana_cputime()
 /*returns a cpu time in seconds */
 {
 int	i;
 i = scalar_scratch(3);
#if NeXT
 sym[i].spec.scalar.f = (float) clock()/( (float) CLK_TCK );
#else
 sym[i].spec.scalar.f = (float) clock()/ 1.e6;
#endif
 return i;
 }
 /*------------------------------------------------------------------------- */
int ana_systime()
 /*returns time since Jan 1, 1970 in dp seconds */
 {
 int	i;
 struct timeval tp;
 struct timezone tzp;
 gettimeofday(&tp, &tzp);
 i = scalar_scratch(4);
 sym[i].spec.scalar.d = (double) tp.tv_sec + .000001* (double) tp.tv_usec;
 return i;
 }
 /*------------------------------------------------------------------------- */
int ana_bang_ctime()
 /*returns current time and date in a string */
 {
 int	i;
 time_t	t;
 i = string_scratch(24);
 t = time(NULL);
 strncpy( (char *) sym[i].spec.array.ptr, ctime( &t ), 24);
 *( (char *) sym[i].spec.array.ptr + 24) = '\0';
 return i;
 }
 /*------------------------------------------------------------------------- */
int ana_ctime(narg,ps)
 /*returns current time and date in a string given a calendar time*/
 int narg,ps[];
 {
 int	result_sym;
 double	sec;
 time_t	t;
 if (double_arg_stat(ps[0], &sec) != 1) return -1;
 t = (time_t) sec;
 result_sym = string_scratch(24);
 strncpy( (char *) sym[result_sym].spec.array.ptr, ctime( &t ), 24);
 *( (char *) sym[result_sym].spec.array.ptr + 24) = '\0';
 return result_sym;
 }
 /*------------------------------------------------------------------------- */
int ana_gmtime(narg,ps)
 /*returns current UT time and date in a string given a calendar time*/
 int narg,ps[];
 {
 int	result_sym;
 double	sec;
 time_t	t;
 struct tm *times;
 if (double_arg_stat(ps[0], &sec) != 1) return -1;
 t = (time_t) sec;
 result_sym = string_scratch(24);
 times = gmtime(&t);
 strncpy( (char *) sym[result_sym].spec.array.ptr, asctime( times ), 24);
 *( (char *) sym[result_sym].spec.array.ptr + 24) = '\0';
 return result_sym;
 }
 /*------------------------------------------------------------------------- */
int ana_mktime(narg,ps)
 /* converts a date to a calendar time using mktime, date is represented
 by a list of integers or a string. The result is returned in double and
 will have any fractional seconds retained. */
 int narg,ps[];
 {
 struct tm times;
 int result_sym;
 float	fsec;
 double	sec;
 if ( sym[ ps[0] ].class == 2 ) {
   char *p, mons[4];
   int	i,j;
   /* decode a string into date parameters */
   p = (char *) sym[ps[0] ].spec.array.ptr;
   sscanf(p, "%s %d,%d %d:%d:%f", mons, &times.tm_mday, &times.tm_year,
    	&times.tm_hour, &times.tm_min, &fsec);
   /* get the month # from the name */
   j = -1;
   for (i=0;i<11;i++)
     //if ( strncmp(tolower(mons), tolower(month_names[i]),3) == 0) j = i;
     if ( strncmp(mons, month_names[i],3) == 0) j = i;
   if (j < 0) {
     printf("bad month in mktime: %s\n", mons);
     return 4;
   }
   times.tm_mon = j;
 } else {
   if (int_arg_stat(ps[0], &times.tm_year) != 1) return -1;
   if (int_arg_stat(ps[1], &times.tm_mon) != 1) return -1;
   times.tm_mon -= 1;
   if (int_arg_stat(ps[2], &times.tm_mday) != 1) return -1;
   if (int_arg_stat(ps[3], &times.tm_hour) != 1) return -1;
   if (int_arg_stat(ps[4], &times.tm_min) != 1) return -1;
   if (float_arg_stat(ps[5], &fsec) != 1) return -1;
 }
 times.tm_isdst = -1;  /* let mktime determine daylight savings */
 times.tm_sec = (time_t) fsec;
 sec = fsec - (float) times.tm_sec;
 if (times.tm_year >= 1900) times.tm_year -= 1900;
 result_sym = scalar_scratch(4);
 sym[result_sym].spec.scalar.d = (double) mktime(&times) + (double) sec;
 return result_sym;
 }
 /*------------------------------------------------------------------------- */
#if defined(SOLARIS) | defined(__sgi) 
time_t my_timegm (struct tm *tm) {
   /* I tried to use the portable suggestion in the timegm man page but because
   setenv is also not available in my problem systems, I am resorting to using
   the timezone variable applied to a mktime result */
   time_t ret, zone;
   extern time_t timezone, altzone, _altzone;
   ret = mktime(tm);
   /* now check for flag */
   if (tm->tm_isdst) zone = _altzone; else zone = timezone;
   ret = mktime(tm) - zone;
   //printf("tm->tm_isdst %d, altzone %d, _altzone %d\n", tm->tm_isdst, timezone, _altzone);
   //printf("mktime result %d, timezone %d, output %d\n", mktime(tm), timezone, ret);
   return ret;
 }
#endif
 /*------------------------------------------------------------------------- */
int ana_timegm(narg,ps)
 /* converts a UTC date to a calendar time using timegm, the date is represented
 by a list of integers or a string. The result is returned in double and
 will have any fractional seconds retained. This is nearly the same as mktime
 and both should be updated and combined. 7/12/2007 */
 int narg,ps[];
 {
 struct tm times;
 int result_sym;
 float	fsec;
 double	sec;
 if ( sym[ ps[0] ].class == 2 ) {
   char *p, mons[4];
   int	i,j;
   /* decode a string into date parameters */
   p = (char *) sym[ps[0] ].spec.array.ptr;
    sscanf(p, "%s %d,%d %d:%d:%f", mons, &times.tm_mday, &times.tm_year,
    	&times.tm_hour, &times.tm_min, &fsec);
   /* get the month # from the name */
   j = -1;
   for (i=0;i<11;i++)
     //if ( strncmp(tolower(mons), tolower(month_names[i]),3) == 0) j = i;
     if ( strncmp(mons, month_names[i],3) == 0) j = i;
   if (j < 0) {
     printf("bad month in mktime: %s\n", mons);
     return 4;
   }
   times.tm_mon = j;
 } else {
   if (int_arg_stat(ps[0], &times.tm_year) != 1) return -1;
   if (int_arg_stat(ps[1], &times.tm_mon) != 1) return -1;
   times.tm_mon -= 1;
   if (int_arg_stat(ps[2], &times.tm_mday) != 1) return -1;
   if (int_arg_stat(ps[3], &times.tm_hour) != 1) return -1;
   if (int_arg_stat(ps[4], &times.tm_min) != 1) return -1;
   if (float_arg_stat(ps[5], &fsec) != 1) return -1;
 }
 times.tm_sec = (time_t) fsec;
 sec = fsec - (float) times.tm_sec;
 if (times.tm_year >= 1900) times.tm_year -= 1900;
 result_sym = scalar_scratch(4);
#if defined(SOLARIS) | defined(__sgi) 
 sym[result_sym].spec.scalar.d = (double) my_timegm(&times) + (double) sec;
#else
 sym[result_sym].spec.scalar.d = (double) timegm(&times) + (double) sec;
#endif
 return result_sym;
 }
 /*------------------------------------------------------------------------- */
int ana_time()
 /*returns current time in a string */
 {
 int	i;
 time_t	t;
 char	*p;
 i = string_scratch(8);
 t = time(NULL);	p = ctime( &t );
 strncpy( (char *) sym[i].spec.array.ptr, p+11, 8);
 *( (char *) sym[i].spec.array.ptr + 8) = '\0';
 return i;
 }
 /*------------------------------------------------------------------------- */
int ana_date()
 /*returns current time in a string */
 {
 int	i;
 time_t	t;
 char	*p, *p2;
 i = string_scratch(12);
 t = time(NULL);	p = ctime( &t );  p2 = (char *) sym[i].spec.array.ptr;
 strncpy( p2, p+4, 6);	strncpy( p2+6, ", ", 2);  strncpy( p2+8, p+20, 4);
 *( (char *) sym[i].spec.array.ptr + 12) = '\0';
 return i;
 }
 /*------------------------------------------------------------------------- */
int ana_symname(narg,ps)
 /*return the symbol name as a string */
 int narg,ps[];
 {
 int	nsym, n;
 char	*s, *s2;
 nsym = ps[0];				/*the target symbol */
 s = find_sym_name(nsym);
 n = strlen(s);
 result_sym = string_scratch(n);
 s2 = (char *) sym[result_sym].spec.array.ptr;
 strcpy(s2, s);
 return result_sym;
 }
 /*------------------------------------------------------------------------- */
int ana_show(narg,ps)
 /*show some info about symbols by number or subname */
 int narg, ps[];
 {
 extern	int	num_internal_syms;
 extern	int	max_sym;
 int	iq, i;
 char	*s, *s2;
 if (narg == 0) return ana_dump(-1, ps);		/* everybody */
 /* if a number passed, show symbol with that number; if a string, find
	 all symbols containing string in name */
 iq = ps[0];
 switch (sym[iq].class) {
 case 8:	iq = class8_to_1(iq);
 case 1:							/*scalar case */
 iq = int_arg(iq);
 if (iq < 0 || iq >= n_symbols) { printf("no such symbol #\n"); return -1; }
 return ana_dump(1, &iq);	/* dump this symbol by number */
 case 2:
						 /* more interesting */
 s = (char *) sym[iq].spec.array.ptr;
 s2 = strsave(s);
 s = s2;
 while (*s2) {*s2 =  toupper( (int) *s2);  s2++; }
 for (i=0;i<max_sym;i++) {
 /* get the name */
 s2 = find_sym_name(i);
 if ( strstr(s2, s) != NULL) { if (ana_dump(1, &i) != 1) break; }
 }
 free(s);
 if (i != max_sym) return -1; else return 1;
 default:
 printf(" only accepts symbol # as an integer or a string for name search\n");
 return -1;
 }
 }
 /*------------------------------------------------------------------------- */
int ana_dump(narg,ps)
 /*show some info about symbols in list */
 int narg, ps[];
 {
 extern	int	num_internal_syms;
 extern	int	max_sym;
 register int	nelem;
 register union types_ptr p1;
 register int	j;
 int i,k,iq,nd,mode,is;
 struct	ahead	*h;
 char	*ptr, *s;
 char *type_string();
 mode = 0;
 if (narg == 0)
	 { mode = 1; narg = max_sym-num_internal_syms; is = num_internal_syms;}
 if (narg <  0)  { mode = 1; narg = max_sym; is = 0;}
 for (i=0;i<narg;i++) {
 if (mode) iq = i+is;  else  iq =ps[i];
 if ( iq < 0 ) return execute_error(54);
 /* get the name */
 s = find_sym_name(iq);
 printf("%4d %10s ",iq,s);
 switch (sym[iq].class)	{		
 case SCALAR:						/*scalar case */
 s = type_string( (int) sym[iq].type);
 printf("scalar, %3s  @0x%08x, value =", s, &sym[iq].spec.scalar.l);
	switch (sym[iq].type) {
	case 0: printf("%10d",sym[iq].spec.scalar.b); break;
	case 1: printf("%10d",sym[iq].spec.scalar.w); break;
	case 2: printf("%10d",sym[iq].spec.scalar.l); break;
	case 3: printf("%14.6e",sym[iq].spec.scalar.f); break;
	case 4: printf("%14.6e",sym[iq].spec.scalar.d); break;
	} printf("\n"); break;			/*end of scalar case */
 case SCAL_PTR:					/*scalar ptr case */
 s = type_string( (int) sym[iq].type);
 printf("scalar ptr, %3s  @0x%08x, value =", s, sym[iq].spec.array.ptr);
	switch (sym[iq].type) {
	case 0: printf("%10d",*(byte *)sym[iq].spec.array.ptr); break;
	case 1: printf("%10d",*(short *)sym[iq].spec.array.ptr); break;
	case 2: printf("%10d",*(int *)sym[iq].spec.array.ptr); break;
	case 3: printf("%14.6e",*(float *)sym[iq].spec.array.ptr); break;
	case 4: printf("%14.6e",*(double *)sym[iq].spec.array.ptr); break;
	} printf("\n"); break;			/*end of scalar ptr case */
 case STRING:							/*string */
	ptr = (char *) sym[iq].spec.array.ptr;
	printf("string, @0x%08x, length: %d value:%.16s\n",ptr,
	sym[iq].spec.array.bstore-1, ptr);
	break;
 case STRING_PTR:						/*strarr */
	h = (struct ahead *) sym[iq].spec.array.ptr;
	nd = h->ndim;
	nelem = 1; for (j=0;j<nd;j++) nelem *= h->dims[j];/*# of elements */
	ptr = (char *) h + sizeof( struct ahead );
	printf("string array @0x%08x, # elem. = %d, (",ptr,nelem);
	for (j=0;j<(nd-1);j++) printf("%d, ",h->dims[j]);
	j=nd-1; printf("%d)\n",h->dims[j]);
	break;
 case ARRAY:							/*array case */
	h = (struct ahead *) sym[iq].spec.array.ptr;
	nd = h->ndim;
	nelem = 1; for (j=0;j<nd;j++) nelem *= h->dims[j];/*# of elements */
	ptr = (char *) h + sizeof( struct ahead );
	s = type_string( (int) sym[iq].type);
	printf("array , %3s  @0x%08x, # elem. = %d, (",s,ptr,nelem);
	for (j=0;j<(nd-1);j++) printf("%d, ",h->dims[j]);
	j=nd-1; printf("%d)\n",h->dims[j]);
	break;
 case SYM_ARR:						/* symarr */
	h = (struct ahead *) sym[iq].spec.array.ptr;
	nd = h->ndim;
	nelem = 1; for (j=0;j<nd;j++) nelem *= h->dims[j];/*# of elements */
	ptr = (char *) h + sizeof( struct ahead );
	printf("symbol array @0x%08x, # elem. = %d, (",ptr,nelem);
	for (j=0;j<(nd-1);j++) printf("%d, ",h->dims[j]);
	j=nd-1; printf("%d)\n",h->dims[j]);
	break;
 case 255:						/* undefined */
	printf("undefined symbol\n"); break;
 default:
	printf("class: %d, type: %d\n", sym[iq].class, sym[iq].type);
	break;
 }						/*end of class switch */
 }						/*end of loop over args */
 return 1;
 }							/*end of ana_dump */
 /*------------------------------------------------------------------------- */
char *type_string(i)
 int	i;
 {
 char	*s;
 switch (i)	{
 case 0: s = "I*1"; break; 
 case 1: s = "I*2"; break;
 case 2: s = "I*4"; break;
 case 3: s = "F*4"; break;
 case 4: s = "F*8"; break;
 default:  execute_error(28); s="???"; break;
 }
 return s;
 }
 /*------------------------------------------------------------------------- */
int ana_zero(narg,ps)
 /*subroutine version,  ZERO, x, [y ...]
 zero symbols in list */
 int narg, ps[];
 {
 int iq, mq, n, type, i;
 union types_ptr	q1;
 char *p;
 for (i=0;i<narg;i++) {
 iq = ps[i];
 /*check if legal to change this symbol */
 if ( iq <= vfix_top ) return execute_error(4);
 /* 6/17/2003 - if undefined, set to a scalar zero, this simplifies
 initializing a bunch of scalars or flags */
 if ( sym[iq].class == 255 || sym[iq].class == 0) {
   int  z=0;
   /* not defined, so do it */
   if (redef_scalar( iq, 2, &z) != 1) return execute_error(13);
 } else
 if (sym[iq].class == 2) {
   /* strings a special case here */
   /* strings are blanked rather than zeroed */
   /* this may not make as much sense in Unix as it seemed to in VMS, may
   want to change */
   mq = sym[iq].spec.array.bstore;  /* note no ahead structure for strings */
   /*mq should now be the # of bytes +1 in the string, get start */
   p = (char *) sym[iq].spec.array.ptr + sizeof( struct ahead );
   mq--;
   while ( mq ) { *p++ = ' ';  mq--; }; *p='\0';	/*but null terminate it */
 } else {
   /* scalars, class 8's, and arrays done here */
   n = get_p_c(&iq, &q1.l);
   if (n <= 0) return -1;
   type = sym[iq].type;
   /* try to zero quickly */
   mq = n * ana_type_size[type];
   bzero( (char *) q1.l, mq);
 }
 }
 return 1;
 }
 /*------------------------------------------------------------------------- */
int ana_zerof(narg,ps)
 /*function version,  x = ZERO(y)
 create array x of same type and size as y and set all elements to zero */
 int narg,ps[];
 {
 int i,iq,jq,mq;
 struct	ahead	*h;
 char *p;
 for (i=0;i<narg;i++) {
 iq =ps[i];
 switch (sym[iq].class) {
 case 8:  iq = class8_to_1(iq);	/*not too sensible, but type is passed on */
 case 1:							/*scalar case */
 /*need an output symbol, are we already a temp ? */
 if ( iq >= temp_base ) result_sym = iq;
   else result_sym = scalar_scratch(sym[iq].type);
	 switch (sym[iq].type) {
	 case 0:sym[result_sym].spec.scalar.b = 0; break;
	 case 1:sym[result_sym].spec.scalar.w = 0; break;
	 case 2:sym[result_sym].spec.scalar.l = 0; break;
	 case 3:sym[result_sym].spec.scalar.f = 0; break;
	 case 4:sym[result_sym].spec.scalar.d = 0; break;
 } return result_sym;				/*end of scalar case */
 case 4:							/*array case */
 /*need an output symbol, are we already a temp ? */
 if ( iq >= temp_base ) result_sym = iq;
   else result_sym = array_clone(iq,sym[iq].type);
 /*try to zero quickly */
 mq = sym[result_sym].spec.array.bstore - sizeof( struct ahead );
 /*mq should now be the # of bytes in the array, get start */
 p = (char *) sym[result_sym].spec.array.ptr + sizeof( struct ahead );
 bzero( p, mq);
 return result_sym;
 case 2:							/*string case */
 /*strings are blanked rather than zeroed */
 /*this may not make as much sense in Unix as it seemed to in VMS, may
 want to change */
 mq = sym[iq].spec.array.bstore;	/*note no ahead structure for strings */
 /*mq should now be the # of bytes +1 in the string, get start */
 p = (char *) sym[iq].spec.array.ptr + sizeof( struct ahead );
 mq--;
 while ( mq ) { *p++ = ' ';  mq--; }; *p='\0';	/*but null terminate it */
 return result_sym;
 }
 }
 return 0;			/*should never get here */
 }
 /*------------------------------------------------------------------------- */
int ana_indgen(narg,ps)
 /*a function, x = indgen(y)
 create array x of same type and size as y and set array elements to their
 index */
 int	narg,ps[];
 {
 register int	j,nelem;
 register union types_ptr p1;
 union	scalar	ct;
 int	iq, mq, nd, id, ind, out, inn, nc, mc;
 char	*ptr;
 struct	ahead	*h;
 iq =ps[0];				/*first argument is symbol */
 switch (sym[iq].class) {
 case 8: iq = class8_to_1(iq);	/*not too sensible, but type is passed on */
 case 1:							/*scalar case */
 /*need an output symbol, are we already a temp ? */
 if ( iq >= temp_base ) result_sym = iq;
  else result_sym = scalar_scratch(sym[iq].type);
	switch (sym[iq].type) {			/*scalars are just zeroed */
	case 0:sym[result_sym].spec.scalar.b = 0; break;
	case 1:sym[result_sym].spec.scalar.w = 0; break;
	case 2:sym[result_sym].spec.scalar.l = 0; break;
	case 3:sym[result_sym].spec.scalar.f = 0; break;
	case 4:sym[result_sym].spec.scalar.d = 0; break;
 }return result_sym;				/*end of scalar case */
 case 4:							/*array case */
 /*need an output symbol, are we already a temp ? */
 if ( iq >= temp_base ) result_sym = iq;
  else result_sym = array_clone(iq,sym[iq].type);
 mq = sym[result_sym].spec.array.bstore - sizeof( struct ahead );
 /*mq should now be the # of bytes in the array, get start */
 h = (struct ahead *) sym[result_sym].spec.array.ptr;
 nd = h->ndim;
 nelem = 1; for (j=0;j<nd;j++) nelem *= h->dims[j];	/*# of elements */
 ptr = (char *) h + sizeof( struct ahead );
 out = 1;
			 /* now find out if there was a second argument */
 if (narg == 1) {
 id = 0;		ind = nelem;	/* fake the general case */
 } else {					/* the 2 arg version */
 id = int_arg( ps[1]);
 if (id >= nd) return execute_error(98);
 ind = h->dims[id];
					 /* get inner and outer counts */
 out = 1;
 if (id < (nd-1) ) for (j=id+1;j<nd;j++) out *=  h->dims[j]; 
 }
 if ( id == 0 ) {			/* no inners, special case */
 switch ( sym[result_sym].type )	{	/*have to switch on type */
 case 2:	p1.l = (int *) ptr;
 while (out--) { ct.l = 0; nc = ind; while (nc--) *p1.l++ = ct.l++;}  break;
 case 3:	p1.f = (float *) ptr;
 while (out--) { ct.f = 0; nc = ind; while (nc--) *p1.f++ = ct.f++;}  break;
 case 0:	p1.b = (byte *) ptr;
 while (out--) { ct.b = 0; nc = ind; while (nc--) *p1.b ++ = ct.b++;}  break;
 case 1:	p1.w = (short *) ptr;
 while (out--) { ct.w = 0; nc = ind; while (nc--) *p1.w++ = ct.w++;}  break;
 case 4:	p1.d = (double *) ptr;
 while (out--) { ct.d = 0; nc = ind; while (nc--) *p1.d++ = ct.d++;}  break;
 }
 } else {					/* need inners as well */
 inn = 1;	for (j=0;j<id;j++) inn *=  h->dims[j];
 switch ( sym[result_sym].type )	{	/*have to switch on type */
 case 2:	p1.l = (int *) ptr;
 while (out--) { ct.l = 0; nc = ind; while (nc--) { mc = inn;
	 while (mc--) {*p1.l++ = ct.l; }  ct.l++;} } break;
 case 3:	p1.f = (float *) ptr;
 while (out--) { ct.f = 0; nc = ind; while (nc--) { mc = inn;
	 while (mc--) {*p1.f++ = ct.f; }  ct.f++;} } break;
 case 0:	p1.b = (byte *) ptr;
 while (out--) { ct.b = 0; nc = ind; while (nc--) { mc = inn;
	 while (mc--) {*p1.b++ = ct.b; }  ct.b++;} } break;
 case 1:	p1.w = (short *) ptr;
 while (out--) { ct.w = 0; nc = ind; while (nc--) { mc = inn;
	 while (mc--) {*p1.w++ = ct.w; }  ct.w++;} } break;
 case 4:	p1.d = (double *) ptr;
 while (out--) { ct.d = 0; nc = ind; while (nc--) { mc = inn;
	 while (mc--) {*p1.d++ = ct.d; }  ct.d++;} } break;
 }
 } return result_sym;
 }						/*end of class switch */
 return execute_error(22);
 }
  /*------------------------------------------------------------------------- */
int ana_complement(narg, ps)
 /* bitwise COMPLEMENT (or NOT) for integers */
 int	narg, ps[];
 {
 union	types_ptr q1,q3;
 int	type, nsym, n;
 nsym = ps[0];
 type = sym[nsym].type;
 /* only for integer types please */
 if (type > 2) return execute_error(29);
 n = get_p_c_r(&nsym, &q1.l, &result_sym, &q3.l);
 if (n <= 0) return -1;
 switch (type) {
  case 1: while (n) { *q3.w++ = ~ (*q1.w++);n--;} break;
  case 2: while (n) { *q3.l++ = ~ (*q1.l++);n--;} break;
  case 0: while (n) { *q3.b++ = ~ (*q1.b++);n--;} break;
 }
 return result_sym;
 }
 /*----------------------------------------------------------------------*/
int ana_not(narg, ps)
 int narg,ps[];
 /* taken from Strous's version, intended to be a logical NOT, the
 bitwise NOT is now called COMPLEMENT (which is a bit long winded
 but I thought COMP might be confusing). Note that following Louie,
 this returns a byte array rather the long used for other logicals;
 maybe they should all return byte arrays ? */
 /* returns a byte 1 for every 0, and 0 for every non-zero element.
   LS 25feb93 */
 {
 int	iq, i, j, type, nd;
 struct ahead	*h;
 register union types_ptr	arg, result;

 iq = ps[0];
 type = sym[iq].type;
 /*if iq is a temp and a byte, we also use it for result, otherwise get a real
 temp */
 switch (sym[iq].class)	{
 case 8: iq = class8_to_1(iq);				/*scalar ptr case */
 case 1:						/*scalar case */
	arg.l = &(sym[iq].spec.scalar.l);
	if ( iq >= temp_base && type == BYTE ) {
	result.b = arg.b; } else {
	iq = scalar_scratch(BYTE);
	result.b = &(sym[iq].spec.scalar.b);
	}
	narg = 1;
	break;
 case 4:
	h = (struct ahead *) sym[iq].spec.array.ptr;
	nd = h->ndim;
	narg = 1; for (j=0;j<nd;j++) narg *= h->dims[j];
	arg.l = (int *) ((char *)h + sizeof(struct ahead));
	/* ck if we can re-use this symbol, if not we get a new iq */
	if ( iq >= temp_base && type == BYTE ) {
	result.b = arg.b; } else {
	iq = array_clone(iq,BYTE);
	h = (struct ahead *) sym[iq].spec.array.ptr;
	result.b = (byte *) ((char *)h + sizeof(struct ahead));
	}
	break;
 default: return execute_error(18); break;
 }

 switch (type)
 { case BYTE:    while (narg--) *result.b++ = (*arg.b++)? 0: 1;  break;
   case WORD:    while (narg--) *result.b++ = (byte)((*arg.w++)? 0: 1);  break;
   case LONG:    while (narg--) *result.b++ = (byte)((*arg.l++)? 0: 1);  break;
   case FLOAT:   while (narg--) *result.b++ = (byte)((*arg.f++)? 0: 1);  break;
   case DOUBLE:  while (narg--) *result.b++ = (byte)((*arg.d++)? 0: 1);  break;
   }
 return iq;
 }
 /*------------------------------------------------------------------------- */
int ana_neg_func(narg,ps)
 /*take the negative of something */
 int narg,ps[];
 {
 union	types_ptr q1,q3;
 int	type, nsym, n;
 nsym = ps[0];
 type = sym[nsym].type;
 n = get_p_c_r(&nsym, &q1.l, &result_sym, &q3.l);
 if (n <= 0) return -1;
 switch (type) {
  case 1: while (n) { *q3.w++ = - (*q1.w++);n--;} break;
  case 2: while (n) { *q3.l++ = - (*q1.l++);n--;} break;
  case 3: while (n) { *q3.f++ = - (*q1.f++);n--;} break;
  case 4: while (n) { *q3.d++ = - (*q1.d++);n--;} break;
  case 0: while (n) { *q3.b++ = (*q1.b++);n--;}
   printf("warning - attempt to take a negative of unsigned byte\n");
   break;
 }
 return result_sym;
 }						/*end of ana_neg_func */
 /*------------------------------------------------------------------------- */
int ana_finite(narg, ps)
 int narg,ps[];
 /* used to check if a value is finite or if all values in an array
 are finite, returns a single I*4 = 1 if everybody is finite and 0 if
 there are NaN's of any type, note that integer types are not checked,
 we just return with a true
 the call is: logical = finite(x) */
 {
 int	iq, type, n;
 union types_ptr	q1;

 num_of_nans = 0;	/* zero this global, it will contain our tally */
 iq = ps[0];
n = get_p_c(&iq, &q1.l);
 if (n <= 0) return -1;
 type = sym[iq].type;
 if (type < 3) return 1;	/* note that this assumes all integers are */
 				/* finite, the 1 is the symbol for a I*4 1 */
 switch (type)
 { case FLOAT: while (n--)  if (finite( *q1.f++) == 0) num_of_nans++;
 	break;
   case DOUBLE: while (n--)  if (finite( *q1.d++) == 0) num_of_nans++;
 	break;
 }
 if (num_of_nans) return 4; else return 1;
 }
 /*------------------------------------------------------------------------- */
int ana_zeronans(narg, ps)
 int narg,ps[];
 /* a subroutine, zero's any NaN's found in the input symbol , also
 counts them */
 {
 int	iq, type, n;
 union types_ptr	q1;

 num_of_nans = 0;	/* zero this global, it will contain our tally */
 iq = ps[0];
 n = get_p_c(&iq, &q1.l);
 if (n <= 0) return -1;
 type = sym[iq].type;
 if (type < 3) return 1;	/* note that this assumes all integers are */
 				/* finite, the 1 is the symbol for a I*4 1 */
 switch (type)
 { case FLOAT: while (n--)  { if (finite(*q1.f) == 0)
 	{ num_of_nans++; *q1.f = 0; } q1.f++; }
 	break;
   case DOUBLE: while (n--) { if (finite( *q1.d) == 0)
   	{ num_of_nans++; *q1.d = 0; } q1.d++; }
 	break;
 }
 return 1;
 }
 /*------------------------------------------------------------------------- */
int ana_isnan(narg, ps)
 int narg,ps[];
 /* the same as finite but returns true if there are any nan's present in
 the variable, so it is the complement of finite
 the call is: logical = isnan(x) */
 {
 /* just call finit and reverse the result */
 if ( ana_finite(narg, ps) == 1 ) return 4; else return 1;
 }
 /*------------------------------------------------------------------------- */
int ana_wherenans(narg,ps)		/* wherenans function */
 /* X=WHERENANS(array), returns locations (1-D index) of NaN's in the
 array or scalar, returns -1 if none found, note that number of NaN's
 is put into !num_of_nans */ 
 int	narg, ps[];
 {
 int	iq, i, j, n, nd, nc, cnt, *p, type, result_sym;
 struct	ahead	*h;
 union	types_ptr q1;
 /* first just get the number of nan's via a call to finite */
 ana_finite(narg, ps);
 
 if (num_of_nans == 0) return 2;  /* if none, we are done */
 nc = get_p_c(ps, &q1.l);
 if (nc <= 0) return -1;
 type = sym[ps[0]].type;
 
 /* we now know that the array had at least one nan and therefore it must
 also be F*4 or F*8, setup a long array to hold the indices */

 n = num_of_nans;
 result_sym = array_scratch(2, 1, &n);		/*for the result */
 h = (struct ahead *) sym[result_sym].spec.array.ptr;
 p = (int *) ((char *)h + sizeof(struct ahead));

 switch (type)
 { case FLOAT:
 	for (i=0;i<nc;i++)  if (finite( *q1.f++) == 0) *p++ = i;
 	break;
   case DOUBLE:
 	for (i=0;i<nc;i++)  if (finite( *q1.d++) == 0) *p++ = i;
 	break;
 }
 return result_sym;
 }
 /*------------------------------------------------------------------------- */
int ana_abs(narg,ps)
 /*take the absolute value of something */
 int narg,ps[];
 {
 int n;
 union	types_ptr q1,q3;
 int	nd, j;
 int	type, iq, nsym;
 nsym = ps[0];
 type = sym[nsym].type;
 n = get_p_c_r(&nsym, &q1.l, &result_sym, &q3.l);
 if (n <= 0) return -1;
 switch (type) {
	case 1: while (n) { if (*q1.w < 0 ) *q3.w++ = - (*q1.w++);
		else *q3.w++ =  (*q1.w++);	n--;} break;
	case 2: while (n) { if (*q1.l < 0 ) *q3.l++ = - (*q1.l++);
		else *q3.l++ =  (*q1.l++);	n--;} break;
	case 3: while (n) { if (*q1.f < 0 ) *q3.f++ = - (*q1.f++);
		else *q3.f++ =  (*q1.f++);	n--;} break;
	case 4: while (n) { if (*q1.d < 0 ) *q3.d++ = - (*q1.d++);
		else *q3.d++ =  (*q1.d++);	n--;} break;
	case 0: while (n) { *q3.b++ =  (*q1.b++);	n--;} break;
 }
 return result_sym;
 }						/* end of ana_abs */
 /*------------------------------------------------------------------------- */
int ana_greg(narg,ps)
 /* compute something for G.M. */
 int narg,ps[];
 {
 /* arguments are b, p, and ph */
 int	n, dim[2];
 float	a, b, p, ph, *pf, x;
 struct	ahead	*h;
 a = float_arg(ps[0]);
 b = float_arg(ps[1]);
 p = float_arg(ps[2]);
 ph = float_arg(ps[3]);
 
 n=500;
 x = .5;
 while (n--)  x = a * (x*(1.-x) + b* sin(6.283185308 * p * x + ph));
 dim[0] = 500;
 result_sym = array_scratch(3, 1, dim);
 h = (struct ahead *) sym[result_sym].spec.array.ptr;
 pf = (float *) ((char *)h + sizeof(struct ahead));
 n = 500;
 while (n--) *pf++ = x = a * (x*(1.-x) + b* sin(6.283185308 * p * x + ph));
 return result_sym;
 }
 /*------------------------------------------------------------------------- */
int ana_total(narg,ps)
 /*total (sum) all elements of an array (also accepts scalars) */
 int narg,ps[];
 {
 /* accumulate in double precision or even modest arrays have bad
 results, several percent error for a 1024x1024 array with equal
 values for example */
 
 float	sum;
 double	dsum;
 int	n;
 union	types_ptr q1;
 int	nd, j;
 int	type;
 int	iq;
 struct	ahead	*h;
 iq = ps[0];
 n = get_p_c(&iq, &q1.l);
 if (n <= 0) return -1;
 type = sym[iq].type;
 sum = 0.;	scount = n;  dsum = 0.0; 
 if (type < 4) result_sym =scalar_scratch(3); else
	{ result_sym =scalar_scratch(4);}
 switch (type) {
	case 1: while (n) { dsum += (double) *q1.w++; n--; } break;
	case 2: while (n) { dsum += (double) *q1.l++; n--; } break;
	case 3: while (n) { dsum += (double) *q1.f++; n--; } break;
	case 4: while (n) { dsum +=          *q1.d++; n--; } break;
	case 0: while (n) { dsum += (double) *q1.b++; n--; } break;
 }
 if (type < 4) sym[result_sym].spec.scalar.f = (float) dsum; else
	sym[result_sym].spec.scalar.d = dsum;
 return result_sym;
 }						/*end of ana_total */
 /*------------------------------------------------------------------------- */
int ana_array_statistics(narg,ps)
 /* returns mean, sdev, and other statistics for an array */
 /* this a subroutine with a variable number of arguments of the
 form:  array_statistics, array, max, min, mean, [sdev, skew, kurtosis]
 where mean, etc are returned values, either float or double */
 int narg,ps[];
 {
 int n, nc, sdev_flag=0, skew_flag=0, kurtosis_flag=0;
 double	fac;
 union	types_ptr q1;
 int	nd, j, *save_ptr;
 int	type;
 int	iq;
 struct	ahead	*h;
 iq = ps[0];
 n = get_p_c(&iq, &q1.l);
 save_ptr = q1.l;
 nc = n;
 if (n <= 0) return -1;
 type = sym[iq].type;
 fac = 1.0/(double) n;
 switch (narg) {
  case 7: kurtosis_flag = 1;
  case 6: skew_flag = 1;
  case 5: sdev_flag = 1;
 }
 /* two major branches, one for doubles and one for everyone else */
 if (type < 4) {
 float	s, sdev, skew, kurtosis, mean, ss, max, min, xq;
 int	imax, imin;
 s=sdev=skew=kurtosis= 0.0;
 /* need a first pass to get the mean value */
 switch (type) {
	case 1:
	  imax = imin = (int) *q1.w;
	  while (n--) { iq = (int) *q1.w++;  s += (float) iq;
	   if (iq > imax) imax = iq; else if (iq < imin) imin = iq; } break;
	case 2:
	  imax = imin =  *q1.l;
	  while (n--) { iq = *q1.l++;  s += (float) iq; 
	   if (iq > imax) imax = iq; else if (iq < imin) imin = iq; } break;
	case 3:
	  max = min =  *q1.f;
	  while (n--) { xq = *q1.f++;  s += xq;
	   if (xq > max) max = xq; else if (xq < min) min = xq; } break;
	case 0:
	  imax = imin = (int)*q1.b;
	  while (n--) { iq = (int) *q1.b++;  s += (float) iq;
	   if (iq > imax) imax = iq; else if (iq < imin) imin = iq; } break;
 }
 mean = s * fac;
 /* finish the mean returned value now */
 redef_scalar(ps[3], 3, &mean);
 if (type == 3) { redef_scalar(ps[1], 3, &max); redef_scalar(ps[2], 3, &min);}
  else { redef_scalar(ps[1], 2, &imax); redef_scalar(ps[2], 2, &imin);}
 /* need a second pass unless only mean was required */
 if (sdev_flag) {
 n = nc;
 q1.l = save_ptr;
 switch (type) {
	case 1: while (n--) {
	 s = (float) *q1.w++ - mean;  ss = s * s;  sdev += ss; 
	 if (skew_flag) { skew += ss*s; if (kurtosis_flag) kurtosis += ss*ss; }
	} break;
	case 2: while (n--) {
	 s = (float) *q1.l++ - mean;  ss = s * s;  sdev += ss; 
	 if (skew_flag) { skew += ss*s; if (kurtosis_flag) kurtosis += ss*ss; }
	} break;
	case 3: while (n--) {
	 s =         *q1.f++ - mean;  ss = s * s;  sdev += ss; 
	 if (skew_flag) { skew += ss*s; if (kurtosis_flag) kurtosis += ss*ss; }
	} break;
	case 0: while (n--) {
	 s = (float) *q1.b++ - mean;  ss = s * s;  sdev += ss; 
	 if (skew_flag) { skew += ss*s; if (kurtosis_flag) kurtosis += ss*ss; }
	} break;
 }
 sdev = sqrt( sdev/(nc-1));
 redef_scalar(ps[4], 3, &sdev);
 if (skew_flag) {	skew = skew/ (nc*pow(sdev, 3));
 			redef_scalar(ps[5], 3, &skew); }
 if (kurtosis_flag) {	kurtosis = kurtosis/ (nc*pow(sdev, 4) - 3.0);
 			redef_scalar(ps[6], 3, &kurtosis); }
 }
 return 1;
 } else {
  /* the doubles case */
 double	s, sdev, skew, kurtosis, mean, ss, max, min, xq;
 s=sdev=skew=kurtosis= 0.0;
 /* need a first pass to get the mean value */
 max = min =  *q1.d;
 while (n--) { xq = *q1.d++;  s += xq;
  if (xq > max) max = xq; else if (xq < min) min = xq; }
 mean = s * fac;
 /* finish the mean returned value now */
 redef_scalar(ps[3], 4, &mean);
 redef_scalar(ps[1], 4, &max); redef_scalar(ps[2], 4, &min);
 /* need a second pass unless only mean was required */
 if (sdev_flag) {
 n = nc;
 q1.l = save_ptr;
 while (n--) {
  s = *q1.d++ - mean;  ss = s * s;  sdev += ss; 
  if (skew_flag) { skew += ss*s; if (kurtosis_flag) kurtosis += ss*ss; }
 }
 
 sdev = sqrt( sdev/(nc-1));
 redef_scalar(ps[4], 4, &sdev);
 if (skew_flag) {	skew = skew/ (nc*pow(sdev, 3));
 			redef_scalar(ps[5], 4, &skew); }
 if (kurtosis_flag) {	kurtosis = kurtosis/ (nc*pow(sdev, 4) - 3.0);
 			redef_scalar(ps[6], 4, &kurtosis); }
 }
 return 1;
 }
 }
 /*------------------------------------------------------------------------- */
int ana_mean(narg,ps)
 /*returns mean of an array, uses total and then divides by scount, which
 was saved by total */
 int	narg,ps[];
 {
 /* 3/2/96, bug fixed which affected means of F*8 arrays */
 int	i, type;
 type = sym[ps[0]].type;
 i = ana_total(narg,ps);
 /*expect i to be a floating scalar containing the sum */
 /* but could also be F*8 */
 switch (type) {
 case 4:  sym[i].spec.scalar.d = sym[i].spec.scalar.d/scount; break;
 default: sym[i].spec.scalar.f = sym[i].spec.scalar.f/scount; break;
 }
 return i;
 }
 /*------------------------------------------------------------------------- */
int ana_dsum(narg,ps)
 /*double total (sum) all elements of an array (also accepts scalars) */
 /*same as ana_total except for a DP accumulator and result */
 int narg,ps[];
 {
 register	double	sum;
 register 	int n;
 union	types_ptr q1;
 int	nd, j;
 int	type;
 int	iq;
 struct	ahead	*h;
 iq =ps[0];
 n = get_p_c(&iq, &q1.l);
 if (n <= 0) return -1;
 type = sym[iq].type;
 sum = 0.;	scount = n;
 result_sym =scalar_scratch(4);
 switch (type) {
	case 1: while (n) { sum += (double) *q1.w++; n--;  } break;
	case 2: while (n) { sum += (double) *q1.l++; n--;  } break;
	case 3: while (n) { sum += (double) *q1.f++; n--;  } break;
	case 4: while (n) { sum += (double) *q1.d++; n--;  } break;
	case 0: while (n) { sum += (double) *q1.b++; n--;  } break;
 }	
 sym[result_sym].spec.scalar.d = sum;
 return result_sym;
 }						/*end of ana_dsum */
#if __APPLE__
 /*------------------------------------------------------------------------- */
int ana_dotproduct(narg,ps)
 /* dot product of 2 arrays */
 int narg,ps[];
 {
 double	sum;
 int n, nn;
 union	types_ptr q1, q2;
 int	nd, j;
 int	type;
 int	iq;
 struct	ahead	*h;
 iq =ps[0];
 nn = get_p_c(&iq, &q1.l);
 if (nn <= 0) return -1;
 type = sym[iq].type;
 if (type != 4) { printf("dotproduct requires F*8\n"); return -1; }
 /* want the second to be of same type and everything as the first */
 iq =ps[1];
 if (type != sym[iq].type) return execute_error(124);
 n = get_p_c(&iq, &q2.l);
 if (n <= 0) return -1;
 if (n != nn) return execute_error(124);

 sum = 0.;
 result_sym =scalar_scratch(4);
 vDSP_dotprD(q1.d, 1, q2.d, 1, &sum, n);
 sym[result_sym].spec.scalar.d = sum;

 return result_sym;
 }						/*end of ana_dsumproduct */
#endif
 /*------------------------------------------------------------------------- */
int ana_dsumproduct(narg,ps)
 /* double total (sum) the product of 2 arrays (dot product) */
 int narg,ps[];
 {
 double	sum;
 int n, nn;
 union	types_ptr q1, q2;
 int	nd, j;
 int	type;
 int	iq;
 struct	ahead	*h;
 iq =ps[0];
 nn = get_p_c(&iq, &q1.l);
 if (nn <= 0) return -1;
 type = sym[iq].type;
 /* want the second to be of same type and everything as the first */
 iq =ps[1];
 if (type != sym[iq].type) return execute_error(124);
 n = get_p_c(&iq, &q2.l);
 if (n <= 0) return -1;
 if (n != nn) return execute_error(124);

 sum = 0.;
 result_sym =scalar_scratch(4);
 switch (type) {
   case 1: while (n) { sum += (double) *q1.w++ * (double) *q2.w++; n--;  } break;
   case 2: while (n) { sum += (double) *q1.l++ * (double) *q2.l++; n--;  } break;
   case 3: while (n) { sum += (double) *q1.f++ * (double) *q2.f++; n--;  } break;
   case 4: while (n) { sum +=  *q1.d++ * *q2.d++; n--;  } break;
   case 0: while (n) { sum += (double) *q1.b++ * (double) *q2.b++; n--;  } break;
 }	
 sym[result_sym].spec.scalar.d = sum;

 return result_sym;
 }						/*end of ana_dsumproduct */
 /*------------------------------------------------------------------------- */
int ana_qlog(narg,ps)
 /* quick and dirty log */
 int narg,ps[];
 {
 double	alogq=2.4786428, clogq=0.7393214;
 double	mantissa, xq;
 int	e, out_type;
 register int	n;
 register union	types_ptr q1, q3;
 union	types_ptr q2;
 int	nd, j;
 int	iq, type;
 struct	ahead	*h;
 iq = ps[0];
 type = sym[iq].type;
 out_type = 3;
 n = get_p_c(&iq, &q2.l);
 q1.f = q2.f;
 if (n <= 0) return -1;

 /*need output symbol, is input available and of correct type ? */
 if (sym[iq].class == 4) {
  if (iq >= temp_base && ((type == out_type)||(out_type == 3 && type ==2)))
	 { result_sym = iq; sym[iq].type = out_type;}
	 else  result_sym = array_clone(iq, out_type);
  h = (struct ahead *) sym[result_sym].spec.array.ptr;
  q3.f = (float *) ((char *)h + sizeof(struct ahead));
 } else {
  result_sym = scalar_scratch(out_type);
  q3.f = &sym[result_sym].spec.scalar.f;
 }

 switch (type) {
  case 0: while (n--) {
    mantissa = frexp((double) *q1.b++, &e);
    *q3.f++ = alogq*(mantissa - 1.0)/(mantissa + clogq) + e;
  }
  break;
   case 1: while (n--) {
    mantissa = frexp((double) *q1.w++, &e);
    *q3.f++ = alogq*(mantissa - 1.0)/(mantissa + clogq) + e;
  }
  break;
  case 2: while (n--) {
    mantissa = frexp((double) *q1.l++, &e);
    *q3.f++ = alogq*(mantissa - 1.0)/(mantissa + clogq) + e;
  }
  break;
  case 3: while (n--) {
    mantissa = frexp((double) *q1.f++, &e);
    *q3.f++ = alogq*(mantissa - 1.0)/(mantissa + clogq) + e;
  }
  break;
  case 4: while (n--) {
    mantissa = frexp( *q1.d++, &e);
    *q3.f++ = alogq*(mantissa - 1.0)/(mantissa + clogq) + e;
  }
  break;
 
 }
 return result_sym;
 }						/*end of ana_qlog */
 /*------------------------------------------------------------------------- */
 /*math functions with 1 argument, all just call math_funcs with appropiate
 code, all have names of form ana_xxx... */
 /*------------------------------------------------------------------------- */
int ana_sin(narg,ps)
 int	narg,ps[];
 {
 return math_funcs(ps[0],0);
 }
 /*------------------------------------------------------------------------- */
int ana_cos(narg,ps)
 int	narg,ps[];
 {
 return math_funcs(ps[0],1);
 }
 /*------------------------------------------------------------------------- */
int ana_tan(narg,ps)
 int	narg,ps[];
 {
 return math_funcs(ps[0],2);
 }
 /*------------------------------------------------------------------------- */
int ana_asin(narg,ps)
 int	narg,ps[];
 {
 return math_funcs(ps[0],3);
 }
 /*------------------------------------------------------------------------- */
int ana_acos(narg,ps)
 int	narg,ps[];
 {
 return math_funcs(ps[0],4);
 }
 /*------------------------------------------------------------------------- */
int ana_atan(narg,ps)
 int	narg,ps[];
 {
 return math_funcs(ps[0],5);
 }
 /*------------------------------------------------------------------------- */
int ana_sinh(narg,ps)
 int	narg,ps[];
 {
 return math_funcs(ps[0],6);
 }
 /*------------------------------------------------------------------------- */
int ana_cosh(narg,ps)
 int	narg,ps[];
 {
 return math_funcs(ps[0],7);
 }
 /*------------------------------------------------------------------------- */
int ana_tanh(narg,ps)
 int	narg,ps[];
 {
 return math_funcs(ps[0],8);
 }
 /*------------------------------------------------------------------------- */
int ana_sqrt(narg,ps)
 int	narg,ps[];
 {
 return math_funcs(ps[0],9);
 }
 /*------------------------------------------------------------------------- */
int ana_cbrt(narg,ps)
 int	narg,ps[];
 {
 return math_funcs(ps[0],10);
 }
 /*------------------------------------------------------------------------- */
int ana_exp(narg,ps)
 int	narg,ps[];
 {
 return math_funcs(ps[0],11);
 }
 /*------------------------------------------------------------------------- */
int ana_expm1(narg,ps)
 int	narg,ps[];
 {
 return math_funcs(ps[0],12);
 }
 /*------------------------------------------------------------------------- */
int ana_log(narg,ps)
 int	narg,ps[];
 {
 return math_funcs(ps[0],13);
 }
 /*------------------------------------------------------------------------- */
int ana_log10(narg,ps)
 int	narg,ps[];
 {
 return math_funcs(ps[0],14);
 }
 /*------------------------------------------------------------------------- */
int ana_log1p(narg,ps)
 int	narg,ps[];
 {
 return math_funcs(ps[0],15);
 }
 /*------------------------------------------------------------------------- */
int ana_erf(narg,ps)
 int	narg,ps[];
 {
 return math_funcs(ps[0],16);
 }
 /*------------------------------------------------------------------------- */
int ana_erfc(narg,ps)
 int	narg,ps[];
 {
 return math_funcs(ps[0],17);
 }
 /*------------------------------------------------------------------------- */
int ana_atan2(narg,ps)
 /*the right way to do atan's, the atan2 function, 2 arguments */
 int	narg,ps[];
 {
 return math_funcs_2f(ps[0],ps[1],18);
 }
 /*------------------------------------------------------------------------- */
int ana_hypot(narg,ps)
 /* takes 2 args like atan2 */
 int	narg,ps[];
 {
 return math_funcs_2f(ps[0],ps[1],26);
 }
 /*------------------------------------------------------------------------- */
int ana_j0(narg,ps)
 int	narg,ps[];
 {
 return math_funcs(ps[0],19);
 }
 /*------------------------------------------------------------------------- */
int ana_j1(narg,ps)
 int	narg,ps[];
 {
 return math_funcs(ps[0],20);
 }
 /*------------------------------------------------------------------------- */
int ana_jn(narg,ps)
 int	narg,ps[];
 {
 return math_funcs_i_f(ps[0],ps[1],25);
 }
 /*------------------------------------------------------------------------- */
int ana_y0(narg,ps)
 int	narg,ps[];
 {
 return math_funcs(ps[0],21);
 }
 /*------------------------------------------------------------------------- */
int ana_y1(narg,ps)
 int	narg,ps[];
 {
 return math_funcs(ps[0],22);
 }
 /*------------------------------------------------------------------------- */
int ana_yn(narg,ps)
 int	narg,ps[];
 {
 return math_funcs_i_f(ps[0],ps[1],24);
 }
 /*------------------------------------------------------------------------- */
int ana_pow(narg,ps)
 int	narg,ps[];
 {
 return math_funcs_2f(ps[0], ps[1],23);
 }
 /*------------------------------------------------------------------------- */
int math_funcs(nsym, code)
 /*general program for floating point functions */
 int	nsym,code;
 {
 register int n;
 register union	types_ptr q3, q1;
 union	types_ptr q2;
 int	j, type, nd, out_type;
 struct	ahead	*h;

 type = sym[nsym].type;
 if ( type == 4 ) out_type = 4; else out_type=3;

 n = get_p_c(&nsym, &q2.l);
 q1.l = q2.l;

 /*need output symbol, is input available and of correct type ? */
 if (sym[nsym].class == 4) {
  if (nsym >= temp_base && ((type == out_type)||(out_type == 3 && type ==2)))
	 { result_sym =nsym; sym[nsym].type=out_type;}
	 else  result_sym = array_clone(nsym,out_type);
  h = (struct ahead *) sym[result_sym].spec.array.ptr;
  q3.f = (float *) ((char *)h + sizeof(struct ahead));
 } else {
  result_sym = scalar_scratch(out_type);
  q3.f = &sym[result_sym].spec.scalar.f;
 }


 /*addresses and count set up, now do the calculations in a loop determined
 by input and output types, only certain combinations possible */
 switch (out_type)	{
 case 3: switch (type) {
  case 3: while (n) { *q3.f++ = (float)(*mc[code])((double) *q1.f++);n--;}
	 break;
  case 0: while (n) { *q3.f++ = (float)(*mc[code])((double) *q1.b++);n--;}
	 break;
  case 1: while (n) { *q3.f++ = (float)(*mc[code])((double) *q1.w++);n--;}
	 break;
  case 2: while (n) { *q3.f++ = (float)(*mc[code])((double) *q1.l++);n--;}
	 break;
 }	break;	/* added break here 10/5/93 */
 case 4:
  while (n) { *q3.d++ = (*mc[code])(*q1.d++);n--;} break;
 }
 return result_sym;
 }						/*end of math_funcs */
 /*------------------------------------------------------------------------- */
int math_funcs_2f(nsym1, nsym2, code)
 /*general program for floating point functions with 2 floating arguments */
 /*messier than the 1 argument case but not as bad as binary ops routines */
 
 /*assumes that the function requires double arguments (most C functions) */
 int	nsym1, nsym2, code;
 {
 register int n;
 register union	types_ptr q1,q2,q3;
 register union	scalar x;
 union scalar y;
 int	j, type1, type2, nd, out_type, nelem;
 struct	ahead	*h;
 type1 = sym[nsym1].type;
 type2 = sym[nsym2].type;
 if ( type1 == 4 || type2 == 4 ) out_type = 4; else out_type=3;
 /*check for class 8's and convert */
 if ( sym[nsym1].class == 8 ) nsym1 = class8_to_1(nsym1);
 if ( sym[nsym2].class == 8 ) nsym2 = class8_to_1(nsym2);
 /*switch on the class combinations*/
 if ( sym[nsym1].class == 1 )
 {							/*first a scalar */
 if ( sym[nsym2].class == 1 )
 {							/*both scalars */
 /*need an output symbol, anybody a temp? */
 /*note that I*4 temps can be used for the F*4 result */
 if (nsym1 >= temp_base)
	{ result_sym =nsym1; sym[nsym1].type=out_type; }
	else  if ( nsym2 >= temp_base)  { result_sym = nsym2;
		sym[result_sym].type = out_type; }
	else result_sym = scalar_scratch(out_type);
 q1.l = &sym[nsym1].spec.scalar.l;
 q2.l = &sym[nsym2].spec.scalar.l;
 q3.l = &sym[result_sym].spec.scalar.l;
 /*because the function call are assumed to require r*8, just convert both
 scalars and do it, return with cast to r*4 if necessary */
 switch (type1) {
	case 4: x.d = *q1.d; break;
	case 3: x.d = (double)*q1.f; break;
	case 2: x.d = (double)*q1.l; break;
	case 1: x.d = (double)*q1.w; break;
	case 0: x.d = (double)*q1.b; break;
 }
 switch (type2) {
	case 4: y.d = *q2.d; break;
	case 3: y.d = (double)*q2.f; break;
	case 2: y.d = (double)*q2.l; break;
	case 1: y.d = (double)*q2.w; break;
	case 0: y.d = (double)*q2.b; break;
 }
 switch (out_type)	{
  case 3: *q3.f = (float)(*mc[code])(x.d,y.d); break;
  case 4: *q3.d = (*mc[code])(x.d,y.d); break;
 }
 return result_sym; }				/*end of both scalars case */
 /*first a scalar but not second, must be an array or illegal */
 if ( sym[nsym2].class == 4 ) {
							/*scalar-array case */
 /*need output symbol, is array available and of correct type ? */
 if ((nsym2>=temp_base) && ((type2==out_type)||(out_type==3 && type2==2)))
	{ result_sym =nsym2; sym[nsym2].type=out_type;}
	else  result_sym = array_clone(nsym2,out_type);
 q1.l = &sym[nsym1].spec.scalar.l;
			/*convert the scalar and put into x, switch on type */
 switch (type1) {
	case 4: x.d = *q1.d; break;
	case 3: x.d = (double)*q1.f; break;
	case 2: x.d = (double)*q1.l; break;
	case 1: x.d = (double)*q1.w; break;
	case 0: x.d = (double)*q1.b; break;
 }
 h = (struct ahead *) sym[nsym2].spec.array.ptr;
 q2.l = (int *) ((char *)h + sizeof(struct ahead));
 nd = h->ndim;
 n = 1; for (j=0;j<nd;j++) n *= h->dims[j];		/*# of elements */
 h = (struct ahead *) sym[result_sym].spec.array.ptr;
 q3.l = (int *) ((char *)h + sizeof(struct ahead));
			/*now do it, switch on out_type and type of array */
 switch (out_type) {
 case 4:	switch (type2) {
  case 4: while (n) {*q3.d++ = (*mc[code])(x.d,*q2.d++); n--;} break;
  case 3: while (n) {*q3.d++ = (*mc[code])(x.d,(double)*q2.f++); n--;} break;
  case 2: while (n) {*q3.d++ = (*mc[code])(x.d,(double)*q2.l++); n--;} break;
  case 1: while (n) {*q3.d++ = (*mc[code])(x.d,(double)*q2.w++); n--;} break;
  case 0: while (n) {*q3.d++ = (*mc[code])(x.d,(double)*q2.b++); n--;} break;
  }  break;	/* added break here 10/5/93 */
 case 3: switch (type2) {
  case 3: while (n) {*q3.f++ = (float)(*mc[code])(x.d,(double)*q2.f++); n--;}
  break;
  case 2: while (n) {*q3.f++ = (float)(*mc[code])(x.d,(double)*q2.l++); n--;}
  break;
  case 1: while (n) {*q3.f++ = (float)(*mc[code])(x.d,(double)*q2.w++); n--;}
  break;
  case 0: while (n) {*q3.f++ = (float)(*mc[code])(x.d,(double)*q2.b++); n--;}
  break;
  }
 }
 return result_sym;
 }
					 /*illegal if here, what was it ? */
 if ( sym[nsym2].class == 2 ) return execute_error(26);
 return execute_error(27);
 }
 if ( sym[nsym1].class == 4 )
 { if ( sym[nsym2].class == 4 ) {			/*both arrays */
 /*need output symbol, anybody available and of correct type ? */
 if ((nsym1>=temp_base) && ((type1==out_type)||(out_type==3 && type1==2)))
	 { result_sym =nsym1; sym[nsym1].type=out_type;}
 else if ((nsym2>=temp_base) && ((type2==out_type)||(out_type==3 && type2==2)))
	 { result_sym =nsym2; sym[nsym2].type=out_type;}
	 else  result_sym = array_clone(nsym1,out_type);
 
 h = (struct ahead *) sym[nsym1].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 */
 nelem = n;
 h = (struct ahead *) sym[nsym2].spec.array.ptr;
 q2.l = (int *) ((char *)h + sizeof(struct ahead));
 nd = h->ndim;
 n = 1; for (j=0;j<nd;j++) n *= h->dims[j];		/*# of elements */
 if ( n != nelem ) return execute_error(24);		/*must match */
 h = (struct ahead *) sym[result_sym].spec.array.ptr;
 q3.l = (int *) ((char *)h + sizeof(struct ahead));
 /*do the casts and calculations */
 switch (out_type)	{
 case 3: switch (type1) {
  case 3: switch (type2) {
   case 3: while (n--) {
    *q3.f++ = (float)(*mc[code])((double) *q1.f++,(double) *q2.f++); }break;
   case 2: while (n--) {
    *q3.f++ = (float)(*mc[code])((double) *q1.f++,(double) *q2.l++); }break;
   case 1: while (n--) {
    *q3.f++ = (float)(*mc[code])((double) *q1.f++,(double) *q2.w++); }break;
   case 0: while (n--) {
    *q3.f++ = (float)(*mc[code])((double) *q1.f++,(double) *q2.b++); }break;
   } break;
  case 2: switch (type2) {
   case 3: while (n--) {
    *q3.f++ = (float)(*mc[code])((double) *q1.l++,(double) *q2.f++); }break;
   case 2: while (n--) {
    *q3.f++ = (float)(*mc[code])((double) *q1.l++,(double) *q2.l++); }break;
   case 1: while (n--) {
    *q3.f++ = (float)(*mc[code])((double) *q1.l++,(double) *q2.w++); }break;
   case 0: while (n--) {
    *q3.f++ = (float)(*mc[code])((double) *q1.l++,(double) *q2.b++); }break;
   } break;
  case 1: switch (type2) {
   case 3: while (n--) {
    *q3.f++ = (float)(*mc[code])((double) *q1.w++,(double) *q2.f++); }break;
   case 2: while (n--) {
    *q3.f++ = (float)(*mc[code])((double) *q1.w++,(double) *q2.l++); }break;
   case 1: while (n--) {
    *q3.f++ = (float)(*mc[code])((double) *q1.w++,(double) *q2.w++); }break;
   case 0: while (n--) {
    *q3.f++ = (float)(*mc[code])((double) *q1.w++,(double) *q2.b++); }break;
   } break;
  case 0: switch (type2) {
   case 3: while (n--) {
    *q3.f++ = (float)(*mc[code])((double) *q1.b++,(double) *q2.f++); }break;
   case 2: while (n--) {
    *q3.f++ = (float)(*mc[code])((double) *q1.b++,(double) *q2.l++); }break;
   case 1: while (n--) {
    *q3.f++ = (float)(*mc[code])((double) *q1.b++,(double) *q2.w++); }break;
   case 0: while (n--) {
    *q3.f++ = (float)(*mc[code])((double) *q1.b++,(double) *q2.b++); }break;
   } break;
  } break;
 case 4: switch (type1) {
  case 4: switch (type2) {
   case 4: while (n--) {
    *q3.d++ = (*mc[code])(*q1.d++,(double) *q2.d++); }break;
   case 3: while (n--) {
    *q3.d++ = (*mc[code])(*q1.d++,(double) *q2.f++); }break;
   case 2: while (n--) {
    *q3.d++ = (*mc[code])(*q1.d++,(double) *q2.l++); }break;
   case 1: while (n--) {
    *q3.d++ = (*mc[code])(*q1.d++,(double) *q2.w++); }break;
   case 0: while (n--) {
    *q3.d++ = (*mc[code])(*q1.d++,(double) *q2.b++); }break;
   } break;
  case 3: /*only type2 = 4 possible */
   while (n--) {
    *q3.d++ = (*mc[code])((double) *q1.f++,*q2.d++); }break;
  case 2:
   while (n--) {
    *q3.d++ = (*mc[code])((double) *q1.l++,*q2.d++); }break;
  case 1:
   while (n--) {
    *q3.d++ = (*mc[code])((double) *q1.w++,*q2.d++); }break;
  case 0:
   while (n--) {
    *q3.d++ = (*mc[code])((double) *q1.b++,*q2.d++); }break;
  } break;
 }
 return result_sym; 				/*end of array array case */
 }
	 /*first an array but not second, must be a scalar or illegal */
 if ( sym[nsym2].class == 1 ) {
							 /*array scalar case */
 /*need output symbol, is array available and of correct type ? */
 if ((nsym1>=temp_base) && ((type1==out_type)||(out_type==3 && type1==2)))
	 { result_sym =nsym1; sym[nsym1].type=out_type;}
	 else  result_sym = array_clone(nsym1,out_type);
			 /*convert the scalar and put into x, switch on type */
 /*note that q1 is used for second (scalar) and q2 for first (array) */
 q1.l = &sym[nsym2].spec.scalar.l;
 switch (type2) {
	 case 4: x.d = *q1.d; break;
	 case 3: x.d = (double)*q1.f; break;
	 case 2: x.d = (double)*q1.l; break;
	 case 1: x.d = (double)*q1.w; break;
	 case 0: x.d = (double)*q1.b; break;
 }
 h = (struct ahead *) sym[nsym1].spec.array.ptr;
 q2.l = (int *) ((char *)h + sizeof(struct ahead));
 nd = h->ndim;
 n = 1; for (j=0;j<nd;j++) n *= h->dims[j];		/*# of elements */
 h = (struct ahead *) sym[result_sym].spec.array.ptr;
 q3.l = (int *) ((char *)h + sizeof(struct ahead));
			 /*now do it, switch on out_type and type of array */
 switch (out_type) {
 case 4:	switch (type1) {
  case 4: while (n) {*q3.d++ = (*mc[code])(*q2.d++,x.d); n--;} break;
  case 3: while (n) {*q3.d++ = (*mc[code])((double)*q2.f++,x.d); n--;} break;
  case 2: while (n) {*q3.d++ = (*mc[code])((double)*q2.l++,x.d); n--;} break;
  case 1: while (n) {*q3.d++ = (*mc[code])((double)*q2.w++,x.d); n--;} break;
  case 0: while (n) {*q3.d++ = (*mc[code])((double)*q2.b++,x.d); n--;} break;
  }	break;		/* added break here 10/5/93 */
 case 3: switch (type1) {
  case 3: while (n) {*q3.f++ = (float)(*mc[code])((double)*q2.f++,x.d); n--;}
  break;
  case 2: while (n) {*q3.f++ = (float)(*mc[code])((double)*q2.l++,x.d); n--;}
  break;
  case 1: while (n) {*q3.f++ = (float)(*mc[code])((double)*q2.w++,x.d); n--;}
  break;
  case 0: while (n) {*q3.f++ = (float)(*mc[code])((double)*q2.b++,x.d); n--;}
  break;
  }
 }
 return result_sym;
 }
					 /*illegal if here, what was it ? */
 if ( sym[nsym2].class == 2 ) return execute_error(26);
 return execute_error(27);
 }
 return -1;					/*shouldn't get here */
 }						/*end of math_funcs_2f */
 /*------------------------------------------------------------------------- */
int math_funcs_i_f(nsym1, nsym2, code)
 /*general program for floating point functions with int and float arguments */
 
 /*assumes that the function requires double arguments (most C functions) */
 int	nsym1, nsym2, code;
 {
 register int n;
 register union	types_ptr q1,q2,q3;
 register union	scalar x;
 union scalar y;
 int	j, type1, type2, nd, out_type, nelem;
 struct	ahead	*h;
 /* always just int the first argument */
 nsym1 = ana_long(1, &nsym1);
 type1 = sym[nsym1].type;
 type2 = sym[nsym2].type;
 out_type = MAX( type2, 3);
 /*check for class 8's and convert */
 if ( sym[nsym1].class == 8 ) nsym1 = class8_to_1(nsym1);
 if ( sym[nsym2].class == 8 ) nsym2 = class8_to_1(nsym2);
 /*switch on the class combinations*/
 if ( sym[nsym1].class == 1 )
 {							/*first a scalar */
 if ( sym[nsym2].class == 1 )
 {							/*both scalars */
 /*need an output symbol, anybody a temp? */
 /*note that I*4 temps can be used for the F*4 result */
 if (nsym1 >= temp_base)
	 { result_sym =nsym1; sym[nsym1].type=out_type; }
	 else  if ( nsym2 >= temp_base)  { result_sym = nsym2;
		 sym[result_sym].type = out_type; }
	 else result_sym = scalar_scratch(out_type);
 x.l  = sym[nsym1].spec.scalar.l;
 q2.l = &sym[nsym2].spec.scalar.l;
 q3.l = &sym[result_sym].spec.scalar.l;
 switch (type2) {
	 case 4: y.d = *q2.d; break;
	 case 3: y.d = (double)*q2.f; break;
	 case 2: y.d = (double)*q2.l; break;
	 case 1: y.d = (double)*q2.w; break;
	 case 0: y.d = (double)*q2.b; break;
 }
 switch (out_type)	{
   case 3: *q3.f = (float)(*mc[code])(x.l,y.d); break;
   case 4: *q3.d = (*mc[code])(x.l,y.d); break;
 }
 return result_sym; }				/*end of both scalars case */
 /*first a scalar but not second, must be an array or illegal */
 if ( sym[nsym2].class == 4 ) {
							 /*scalar-array case */
 /*need output symbol, is array available and of correct type ? */
 if ((nsym2>=temp_base) && ((type2==out_type)||(out_type==3 && type2==2)))
	 { result_sym =nsym2; sym[nsym2].type=out_type;}
	 else  result_sym = array_clone(nsym2,out_type);
 x.l = sym[nsym1].spec.scalar.l;
 h = (struct ahead *) sym[nsym2].spec.array.ptr;
 q2.l = (int *) ((char *)h + sizeof(struct ahead));
 nd = h->ndim;
 n = 1; for (j=0;j<nd;j++) n *= h->dims[j];		/*# of elements */
 h = (struct ahead *) sym[result_sym].spec.array.ptr;
 q3.l = (int *) ((char *)h + sizeof(struct ahead));
			 /*now do it, switch on out_type and type of array */
 switch (out_type) {
 case 4:	switch (type2) {
  case 4: while (n) {*q3.d++ = (*mc[code])(x.l,*q2.d++); n--;} break;
  case 3: while (n) {*q3.d++ = (*mc[code])(x.l,(double)*q2.f++); n--;} break;
  case 2: while (n) {*q3.d++ = (*mc[code])(x.l,(double)*q2.l++); n--;} break;
  case 1: while (n) {*q3.d++ = (*mc[code])(x.l,(double)*q2.w++); n--;} break;
  case 0: while (n) {*q3.d++ = (*mc[code])(x.l,(double)*q2.b++); n--;} break;
  }	break;		/* added break here 10/5/93 */
 case 3: switch (type2) {
  case 3: while (n) {*q3.f++ = (float)(*mc[code])(x.l,(double)*q2.f++); n--;}
  break;
  case 2: while (n) {*q3.f++ = (float)(*mc[code])(x.l,(double)*q2.l++); n--;}
  break;
  case 1: while (n) {*q3.f++ = (float)(*mc[code])(x.l,(double)*q2.w++); n--;}
  break;
  case 0: while (n) {*q3.f++ = (float)(*mc[code])(x.l,(double)*q2.b++); n--;}
  break;
  }
 }
 return result_sym;
 }
					 /*illegal if here, what was it ? */
 if ( sym[nsym2].class == 2 ) return execute_error(26);
 return execute_error(27);
 }
 if ( sym[nsym1].class == 4 )
 { if ( sym[nsym2].class == 4 ) {			/*both arrays */
 /*need output symbol, anybody available and of correct type ? */
 if ((nsym1>=temp_base) && ((type1==out_type)||(out_type==3 && type1==2)))
	 { result_sym =nsym1; sym[nsym1].type=out_type;}
 else if ((nsym2>=temp_base) && ((type2==out_type)||(out_type==3 && type2==2)))
	 { result_sym =nsym2; sym[nsym2].type=out_type;}
	 else  result_sym = array_clone(nsym1,out_type);
 
 h = (struct ahead *) sym[nsym1].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 */
 nelem = n;
 h = (struct ahead *) sym[nsym2].spec.array.ptr;
 q2.l = (int *) ((char *)h + sizeof(struct ahead));
 nd = h->ndim;
 n = 1; for (j=0;j<nd;j++) n *= h->dims[j];		/*# of elements */
 if ( n != nelem ) return execute_error(24);		/*must match */
 h = (struct ahead *) sym[result_sym].spec.array.ptr;
 q3.l = (int *) ((char *)h + sizeof(struct ahead));
 /*do the casts and calculations */
 switch (out_type)	{
 case 3:
  switch (type2) {
   case 3: while (n--) {
    *q3.f++ = (float)(*mc[code])((double) *q1.l++,(double) *q2.f++); }break;
   case 2: while (n--) {
    *q3.f++ = (float)(*mc[code])((double) *q1.l++,(double) *q2.l++); }break;
   case 1: while (n--) {
    *q3.f++ = (float)(*mc[code])((double) *q1.l++,(double) *q2.w++); }break;
   case 0: while (n--) {
    *q3.f++ = (float)(*mc[code])((double) *q1.l++,(double) *q2.b++); }break;
   } break;
 case 4:
  switch (type2) {
   case 4: while (n--) {
    *q3.d++ = (*mc[code])(*q1.l++,(double) *q2.d++); }break;
   case 3: while (n--) {
    *q3.d++ = (*mc[code])(*q1.l++,(double) *q2.f++); }break;
   case 2: while (n--) {
    *q3.d++ = (*mc[code])(*q1.l++,(double) *q2.l++); }break;
   case 1: while (n--) {
    *q3.d++ = (*mc[code])(*q1.l++,(double) *q2.w++); }break;
   case 0: while (n--) {
    *q3.d++ = (*mc[code])(*q1.l++,(double) *q2.b++); }break;
   } break;
 }
 return result_sym; 				/*end of array array case */
 }
	 /*first an array but not second, must be a scalar or illegal */
 if ( sym[nsym2].class == 1 ) {
							 /*array scalar case */
 /*need output symbol, is array available and of correct type ? */
 if ((nsym1>=temp_base) && ((type1==out_type)||(out_type==3 && type1==2)))
	 { result_sym =nsym1; sym[nsym1].type=out_type;}
	 else  result_sym = array_clone(nsym1,out_type);
			 /*convert the scalar and put into x, switch on type */
 /*note that q1 is used for second (scalar) and q2 for first (array) */
 q1.l = &sym[nsym2].spec.scalar.l;
 switch (type2) {
	 case 4: x.d = *q1.d; break;
	 case 3: x.d = (double)*q1.f; break;
	 case 2: x.d = (double)*q1.l; break;
	 case 1: x.d = (double)*q1.w; break;
	 case 0: x.d = (double)*q1.b; break;
 }
 h = (struct ahead *) sym[nsym1].spec.array.ptr;
 q2.l = (int *) ((char *)h + sizeof(struct ahead));
 nd = h->ndim;
 n = 1; for (j=0;j<nd;j++) n *= h->dims[j];		/*# of elements */
 h = (struct ahead *) sym[result_sym].spec.array.ptr;
 q3.l = (int *) ((char *)h + sizeof(struct ahead));
			 /*now do it, switch on out_type and type of array */
 switch (out_type) {
 case 4:  while (n) {*q3.d++ = (*mc[code])((double)*q2.l++,x.d); n--;} break;
 case 3:  while (n) {*q3.f++ = (float)(*mc[code])((double)*q2.l++,x.d); n--;}
  break;
 }
 return result_sym;
 }
					 /*illegal if here, what was it ? */
 if ( sym[nsym2].class == 2 ) return execute_error(26);
 return execute_error(27);
 }
 return -1;					/*shouldn't get here */
 }						/*end of math_funcs_i_f */
 /*------------------------------------------------------------------------- */
int ana_therm2temp(narg,ps)
 int	narg,ps[];
 {
 return spec_funcs(ps[0],0);
 }
 /*------------------------------------------------------------------------- */
int ana_patan(narg,ps)
 int	narg,ps[];
 {
 return spec_funcs(ps[0],1);
 }
 /*------------------------------------------------------------------------- */
int spec_funcs(nsym, code)
 /*general program for special integer functions */
 int	nsym,code;
 {
 register int n;
 register union	types_ptr q3, q1;
 union	types_ptr q2;
 int	j, type, nd, out_type;
 int	(*ifun)() = mcint[code];
 struct	ahead	*h;

 type = sym[nsym].type;
 out_type=2;

 n = get_p_c(&nsym, &q2.l);
 q1.l = q2.l;

 /*need output symbol, is input available and of correct type ? */
 if (sym[nsym].class == 4) {
  if ( nsym >= temp_base && (type == out_type) )
	 { result_sym =nsym; }
	 else  result_sym = array_clone(nsym, out_type);
  h = (struct ahead *) sym[result_sym].spec.array.ptr;
  q3.l = (int *) ((char *)h + sizeof(struct ahead));
 } else {
  result_sym = scalar_scratch(out_type);
  q3.l = &sym[result_sym].spec.scalar.l;
 }


 /*addresses and count set up, now do the calculations in a loop,
 out_type is always 2 here */
 switch (type) {
  case 4: while (n) { *q3.l++ = (int)(*ifun)((int) *q1.d++);n--;}
	 break;
  case 3: while (n) { *q3.l++ = (int)(*ifun)((int) *q1.f++);n--;}
	 break;
  case 0: while (n) { *q3.l++ = (int)(*ifun)((int) *q1.b++);n--;}
	 break;
  case 1: while (n) { *q3.l++ = (int)(*ifun)((int) *q1.w++);n--;}
	 break;
  case 2: while (n) { *q3.l++ = (int)(*ifun)((int) *q1.l++);n--;}
	 break;
 }
 return result_sym;
 }						/*end of math_funcs */
 /*------------------------------------------------------------------------- */