/* internal ana subroutines and functions */ /*this is file fun2.c */ #include #include #include #include #include #include #include "ana_structures.h" #include "defs.h" /* for SGI only (?) */ #if __sgi #include #endif #define string_func_setup \ if ( sym[ ps[0] ].class != STRING ) return execute_error(70);\ if ( sym[ ps[1] ].class != STRING ) return execute_error(70);\ p1 = (char *) sym[ps[0] ].spec.array.ptr;\ p2 = (char *) sym[ps[1] ].spec.array.ptr; #define linesize 256 extern struct sym_desc sym[]; extern struct sym_list *subr_sym_list[]; extern int temp_base; extern int ana_float(); extern char *find_sym_name(int iq); extern int ana_type_size[]; extern int vfix_top, num_ana_subr, next_user_subr_num; extern int num_ana_func, next_user_func_num; extern float float_arg(); extern double double_arg(); extern int get_p_c(), get_p_c_r(); extern char *strsave(); union types_ptr { byte *b; short *w; int *l; float *f; double *d;}; static int result_sym , type, detrend_flag; static char line[linesize]; /* a temp for output and istring */ float gsmooth_width = 2.3; /* limit for gsmooth kernel, good for .005, use 4 for 1.e-7 */ /*------------------------------------------------------------------------- */ int ana_runsum(narg,ps) /*generate a running sum */ /*this version floats the array first, it might be faster to use a switch so that the conversion to float is done in the inner loop, also support double*/ int narg, ps[]; { float sum; int n; union types_ptr q1,q2; int nd, j; int iq; struct ahead *h; iq = ps[0]; /*float the input */ iq = ana_float(narg,ps); n = get_p_c_r(&iq, &q1, &result_sym, &q2); if (n <= 0) return -1; sum = 0.; while (n) { sum += *q1.f++; *q2.f++ = sum; n--; } return result_sym; } /*end of ana_runsum */ /*------------------------------------------------------------------------- */ int ana_sdev(narg,ps) /*compute the sample standard deviation of an array */ /*this version floats the array first, it might be faster to use a switch so that the conversion to float is done in the inner loop, also support double*/ int narg,ps[]; { float sum, sq, xq; int n; union types_ptr q1; int nd, j, nsave; int iq; float *psave; iq = ps[0]; /*float the input */ iq = ana_float(narg,ps); n = get_p_c(&iq, &q1); if (n <= 0) return -1; result_sym =scalar_scratch(3); /*first pass through to get the mean or the roundoffs will get us */ sum = 0.; nsave = n; psave = q1.f; while (n) { sum += *q1.f++; n--; } sum = sum/( (float) nsave); sq = 0.0; n = nsave; q1.f = psave; while (n) { xq = (*q1.f++) - sum; sq += (xq * xq); n--; } /*mean in sum and sum of squares in sq */ if (nsave > 1 ) xq = sqrt( sq / ((float) (nsave -1) )); else xq = 0.0; sym[result_sym].spec.scalar.f = xq; return result_sym; } /*end of ana_sdev */ /*------------------------------------------------------------------------- */ int ana_swab(narg,ps) /*swap bytes in an array or even an scalar */ /*this is a subroutine and swaps the input symbol (not a copy) */ int narg,ps[]; { int iq, n, nd, j; union types_ptr q1; struct ahead *h; iq = ps[0]; /*just 1 argument */ /*check if legal to change this symbol */ if ( iq <= vfix_top ) return execute_error(4); type = sym[iq].type; n = get_p_c(&iq, &q1); if (n <= 0) return -1; n = n * ana_type_size[type]; /* get byte count */ swapb( (char *) q1.l, n); return 1; } /*------------------------------------------------------------------------- */ int ana_swapl(narg,ps) /*reverse bytes in each I*4 part of an array or scalar */ /*this is a subroutine and swapl's the input symbol (not a copy) */ int narg,ps[]; { int iq, n, nd, j; union types_ptr q1; struct ahead *h; iq = ps[0]; /*just 1 argument */ /*check if legal to change this symbol */ if ( iq <= vfix_top ) return execute_error(4); type = sym[iq].type; n = get_p_c(&iq, &q1); if (n <= 0) return -1; n = n * ana_type_size[type]; /* get byte count */ swapl( (char *) q1.l, n/4); return 1; } /*------------------------------------------------------------------------- */ int ana_endian(narg,ps) /*reverse byte order depending on type */ /*this is a subroutine and swap's the input symbol (not a copy) */ int narg,ps[]; { int iq, n, nd, j; char *q; struct ahead *h; iq = ps[0]; /*just 1 argument */ /*check if legal to change this symbol */ if ( iq <= vfix_top ) return execute_error(4); type = sym[iq].type; n = get_p_c(&iq, &q); if (n <= 0) return -1; n = n * ana_type_size[type]; /* get byte count */ if (type == 1) swapb(q, n); if (type == 2 || type ==3) swapl(q, n/4); if (type == 4) swapd(q, n/8); return 1; } /*------------------------------------------------------------------------- */ int ana_wait(narg,ps) /* wait the specified # of seconds (floating point) */ int narg,ps[]; { float xq; struct timeval tval; xq = float_arg(ps[0]); xq = MAX(0.0, xq); tval.tv_sec = (unsigned int) xq; xq = (xq - tval.tv_sec) * 1.e6; tval.tv_usec = (unsigned int) xq; select(0, NULL, NULL, NULL,&tval); return 1; } /*------------------------------------------------------------------------- */ void wait_sec(xq) /* for internal use */ float xq; { struct timeval tval; xq = MAX(0.0, xq); tval.tv_sec = (unsigned int) xq; xq = (xq - tval.tv_sec) * 1.e6; tval.tv_usec = (unsigned int) xq; select(0, NULL, NULL, NULL,&tval); } /*------------------------------------------------------------------------- */ int ana_smooth(narg,ps) /* returns a smooth version of input, boxcar smooth */ /* ps[0] is input array and ps[1] is smooth width (always rounded to odd) */ int narg,ps[]; { float sum, *pt1, *pt2, wq; int n; union types_ptr q1,q2; int nd, j, width; int iq, oldtype, nxq, outer; struct ahead *h; iq = ps[0]; if ( sym[iq].class != 4 ) return execute_error(66); oldtype = sym[iq].type; /*original type */ /*float the input */ iq = ana_float(1,ps); if (int_arg_stat(ps[1], &width) != 1) return -1; if ( width <= 1 ) return ps[0]; /* just return original if width too small*/ if (width%2 == 0 ) width += 1; /*make next odd value if even */ h = (struct ahead *) sym[iq].spec.array.ptr; q1.l = (int *) ((char *)h + sizeof(struct ahead)); nd = h->ndim; nxq = h->dims[0]; outer = 1; if (nd > 1) for(j=1;jdims[j]; result_sym = array_clone(iq,3); h = (struct ahead *) sym[result_sym].spec.array.ptr; q2.l = (int *) ((char *)h + sizeof(struct ahead)); wq = 1.0/width; /* we do 1-D smooth over inner dimension, looping over any outer dimensions*/ if ( width >= nxq ) { /* a nearly trivial case */ while (outer--) { n = nxq; sum = 0.0; while (n--) sum += *q1.f++; sum = sum / nxq; n = nxq; while (n--) *q2.f++ = sum; } } else { /*normal case */ while (outer--) { sum = 0.0; n = width; pt1 = q1.f; while (n--) sum += *q1.f++; /* starting sum */ sum = sum * wq; n = width/2; while (n--) *q2.f++ = sum; n = nxq - width; /*middle section */ while (n--) { *q2.f++ = sum; sum += wq * ( *q1.f++ - *pt1++ ); } n = width/2 + 1; /*end section */ while (n--) { *q2.f++ = sum; } } } /* convert back to oldtype if it wasn't floating */ if ( oldtype == 3 ) return result_sym; switch (oldtype) { case 0: return ana_byte(1,&result_sym); case 1: return ana_word(1,&result_sym); case 2: return ana_long(1,&result_sym); case 4: return ana_double(1,&result_sym); } return execute_error(3); } /*------------------------------------------------------------------------- */ int ana_bsmooth(narg,ps) /* binomial filter */ /* ps[0] is input array and ps[1] is smooth width (always rounded to odd) */ int narg,ps[]; { int n, nmo, npass, nps; union types_ptr q1, q2; int nd, j, width; int iq, type, nxq, outer; struct ahead *h; iq = ps[0]; if ( sym[iq].class != 4 ) return execute_error(66); type = sym[iq].type; /*original type */ if (int_arg_stat(ps[1], &width) != 1) return -1; if ( width <= 1 ) return ps[0]; /* just return original if width too small*/ if (width%2 == 0 ) width += 1; /*make next odd value if even */ npass = (width - 1)/2; h = (struct ahead *) sym[iq].spec.array.ptr; q1.l = (int *) ((char *)h + sizeof(struct ahead)); nd = h->ndim; nxq = h->dims[0]; outer = 1; if (nd > 1) for(j=1;jdims[j]; result_sym = array_clone(iq, type); h = (struct ahead *) sym[result_sym].spec.array.ptr; q2.l = (int *) ((char *)h + sizeof(struct ahead)); /* we do 1-D smooth over inner dimension, looping over any outer dimensions*/ nmo = nxq - 1; switch (type) { case 0: { int a1, a2, t1, t2, acc; byte *p = q1.b, *q = q2.b; while (outer--) { nps = npass; while (nps-- > 0) { a1 = (int) *p++; a2 = 2*a1; n = nmo; while (n--) { acc = *p++; t1 = acc; acc += a1; a1 = t1; t2 = acc; acc += a2; a2 = t2; *q++ = (byte) (acc>>2); } *q++ = (byte) ((a1 + a1 + a2)>>2); /* for additional passes, we smooth over previous result */ p = q = q2.b; } q1.b += nxq; q2.b += nxq; p = q1.b; q = q2.b; } } break; case 1: { int a1, a2, t1, t2, acc; short *p = q1.w, *q = q2.w; while (outer--) { nps = npass; while (nps-- > 0) { a1 = (int) *p++; a2 = 2*a1; n = nmo; while (n--) { acc = *p++; t1 = acc; acc += a1; a1 = t1; t2 = acc; acc += a2; a2 = t2; *q++ = (short) (acc>>2); } *q++ = (short) ((a1 + a1 + a2)>>2); /* for additional passes, we smooth over previous result */ p = q = q2.w; } q1.w += nxq; q2.w += nxq; p = q1.w; q = q2.w; } } break; case 2: /* 32 bit image, need to use a larger accumulator or suffer limitations, long long type used, probably be slow on most machines */ { long long a1, a2, t1, t2, acc; int *p = q1.l, *q = q2.l; while (outer--) { nps = npass; while (nps-- > 0) { a1 = (int) *p++; a2 = 2*a1; n = nmo; while (n--) { acc = *p++; t1 = acc; acc += a1; a1 = t1; t2 = acc; acc += a2; a2 = t2; *q++ = (short) (acc>>2); } *q++ = (short) ((a1 + a1 + a2)>>2); /* for additional passes, we smooth over previous result */ p = q = q2.l; } q1.l += nxq; q2.l += nxq; p = q1.l; q = q2.l; } } break; case 3: /* printf("float case, nmo = %d, outer = %d\n", nmo, outer); */ { float a1, a2, t1, t2, acc; float *p = q1.f, *q = q2.f; while (outer--) { nps = npass; while (nps-- > 0) { a1 = *p++; a2 = 2*a1; n = nmo; while (n--) { acc = *p++; t1 = acc; acc += a1; a1 = t1; t2 = acc; acc += a2; a2 = t2; *q++ = acc * 0.25; } *q++ = (a1 + a1 + a2) * 0.25; /* for additional passes, we smooth over previous result */ p = q = q2.f; } q1.f += nxq; q2.f += nxq; p = q1.f; q = q2.f; } } break; case 4: { double a1, a2, t1, t2, acc; double *p = q1.d, *q = q2.d; while (outer--) { nps = npass; while (nps-- > 0) { a1 = *p++; a2 = 2*a1; n = nmo; while (n--) { acc = *p++; t1 = acc; acc += a1; a1 = t1; t2 = acc; acc += a2; a2 = t2; *q++ = acc * 0.25; } *q++ = (a1 + a1 + a2) * 0.25; /* for additional passes, we smooth over previous result */ p = q = q2.d; } q1.d += nxq; q2.d += nxq; p = q1.d; q = q2.d; } } break; } return result_sym; } /*------------------------------------------------------------------------- */ int ana_gsmooth(int narg,int ps[]) /* smooth input array with a gaussian */ /* ps[0] is input array and ps[1] is fwhm width */ /* add an optional mask that limits the smooth to the mask */ { extern int scrat[NSCRAT]; float sum, *pt1, *pt2, *pt3, wq, sumg, xq; int n; union types_ptr q1, q2; int nd, j, n2, ng, i2; int iq, oldtype, nx, ny, outer; float width, *gkern, gsum; struct ahead *h; iq = ps[0]; if ( sym[iq].class != 4 ) return execute_error(66); /* 7/28/2010 - don't convert back to original type, always output float, this reduces conversion overhead when used in gsmooth2 or followed by gsmoothy */ //oldtype = sym[iq].type; /*original type */ /* still float the input, may want to only float pixels used for mask case (someday) */ iq = ana_float(1,&iq); width = 0.6005612 * float_arg( ps[1] ); /* convert fwhm to gaussian*/ h = (struct ahead *) sym[iq].spec.array.ptr; q1.l = (int *) ((char *)h + sizeof(struct ahead)); nd = h->ndim; nx = h->dims[0]; outer = 1; if (nd > 1) { ny = h->dims[1]; for(j=1;jdims[j]; } result_sym = array_clone(iq,3); h = (struct ahead *) sym[result_sym].spec.array.ptr; q2.l = (int *) ((char *)h + sizeof(struct ahead)); /* check for trivial cases, just return original */ if ( width < 0.25 ) return ps[0]; /* set up the gaussian convolution kernel */ n2 = MIN( (int) (gsmooth_width * width), nx); n2 = MAX( n2, 1); ng = 2 * n2 +1; if ( (ng) <= NSCRAT) gkern = (float *) scrat; else gkern = (float *) malloc( (ng) * sizeof(float)); n = ng; wq = (float) (-n2); pt3 = gkern; gsum = 0.0; while (n--) { sum = wq/width; wq++; xq = exp(-sum*sum); *pt3++ = xq; gsum += xq; } /* now check if we have a mask */ if (narg <= 2) { /* no mask case, we do 1-D smooth over inner dimension, looping over any outer dimensions*/ // /*float the input */ // iq = ps[0]; // iq = ana_float(1,&iq); // h = (struct ahead *) sym[iq].spec.array.ptr; // q1.f = (float *) ((char *)h + sizeof(struct ahead)); pt2 = q2.f; while (outer--) { /* zone 1 */ n = MIN(n2, nx); for (j=0;jndim; if (nd == 1) { /* can't be mask, check if legal index array */ int ntotal, *indptr; ntotal = get_1d_array(iq, LONG, &indptr); printf("GSMOOTH - index array version not supported yet\n"); return -1; } else { /* mask possibility, check */ int ntotal, n, m, nc1, nc2, nc3; byte *maskbase, *maskptr, *maskptr2; printf("mask case in GSMOOTH\n"); ntotal = get_2d_array(iq, BYTE, (int **) &maskbase, &n, &m); if (ntotal <= 0 || n != nx || m != ny) { printf("GSMOOTH - mask array must be 2-D byte array and dimensions must match image\n"); return -1; } /* OK, a mask it is then */ maskptr = maskptr2 = maskbase; bzero(q2.f, ntotal * sizeof(float)); pt2 = q2.f; nc1 = MIN(n2, nx); nc2 = MAX(0, nx - 2 * n2); //nc3 = n2; if ( nx < (2 * n2) ) nc3 = MAX(0, nx - n2); while (outer--) { /* zone 1 */ n = nc1; for (j=0;jndim; if (nd != 2) { fprintf(stderr, "GSMOOTHY - array must be 2-D, not %d\n", nd); return -1; } ny = h->dims[1]; nx = inner = h->dims[0]; result_sym = array_clone(iq,3); h = (struct ahead *) sym[result_sym].spec.array.ptr; q2.l = (int *) ((char *)h + sizeof(struct ahead)); /* check for trivial cases, just return original */ if ( width < 0.25 ) return ps[0]; /* set up the gaussian convolution kernel */ n2 = MIN( (int) (gsmooth_width * width), ny); n2 = MAX( n2, 1); ng = 2 * n2 +1; if ( (ng) <= NSCRAT) gkern = (float *) scrat; else gkern = (float *) malloc( (ng)* sizeof(float)); n = ng; wq = (float) (-n2); pt3 = gkern; gsum = 0.0; while (n--) { sum = wq/width; wq++; xq = exp(-sum*sum); *pt3++ = xq; gsum += xq; } /* now check if we have a mask */ if (narg <= 2) { /* we do 1-D smooth over outer dimension, looping over any inner dimensions*/ /* row buffers are setup or we get badly punished on memory access */ rowbuf = (float *) malloc( (ny)* sizeof(float)); outrowbuf = (float *) malloc( (ny)* sizeof(float)); while (inner--) { pr = rowbuf; pt1 = q1.f; n = ny; while (n--) { *pr++ = *pt1; pt1 += nx; } pt1 = rowbuf; /* now pt1 is used the same as in gsmooth */ pt2 = outrowbuf; /* zone 1 */ n = MIN(n2, ny); for (j=0;jndim; if (nd == 1) { /* can't be mask, check if legal index array */ int ntotal, *indptr; ntotal = get_1d_array(iq, LONG, &indptr); printf("GSMOOTHY - index array version not supported yet\n"); return -1; } else { /* mask possibility, check */ int ntotal, n, m, nc1, nc2, nc3; byte *maskbase, *maskptr, *maskptr2; printf("mask case in GSMOOTH\n"); ntotal = get_2d_array(iq, BYTE, (int **) &maskbase, &n, &m); if (ntotal <= 0 || n != nx || m != ny) { printf("GSMOOTH - mask array must be 2-D byte array and dimensions must match image\n"); return -1; } /* OK, a mask it is then */ maskptr = maskptr2 = maskbase; bzero(q2.f, ntotal * sizeof(float)); pt2 = q2.f; nc1 = MIN(n2, ny); nc2 = MAX(0, ny - 2 * n2); nc3 = n2; if ( ny < (2 * n2) ) nc3 = MAX(0, ny - n2); while (inner--) { /* zone 1 */ n = nc1; for (j=0;jndim; nxq = h->dims[0]; outerx = 1; if (nd > 1) for(j=1;jdims[j]; /* now y */ iq = ps[1]; if ( sym[iq].class != 4 ) return execute_error(66); iq = ana_float(1,&iq); h = (struct ahead *) sym[iq].spec.array.ptr; y = (float *) ((char *)h + sizeof(struct ahead)); nd = h->ndim; nyq = h->dims[0]; outery = 1; dim[0] = nyq; if (nd > 1) for(j=1;jdims[j]; dim[j] = h->dims[j]; } /* nxq must = nyq and outerx must = outery or 1 */ if ( nxq != nyq || ( outerx != outery && outerx != 1) ) return execute_error(103); /* now f, make a match for y array */ iq = ps[2]; if (redef_array(iq, 3, nd, dim) != 1) return -1; h = (struct ahead *) sym[iq].spec.array.ptr; f = (float *) ((char *)h + sizeof(struct ahead)); /* check for presence of fixed parameters */ if ( (gfit_fix_a && narg < 4) || (gfit_fix_b && narg < 5) || (gfit_fix_x0 && narg < 6) || (gfit_fix_del && narg < 7) ) { printf("GFIT: missing fixed parameter\n"); return -1; } /* a place to point to for optional args */ pa = &guess[0]; pb = &guess[1]; px0 = &guess[2]; pdel = &guess[3]; inc_a = inc_b = inc_x0 = inc_del = 0; update_a = update_b = update_x0 = update_del = 0; /* the next 4 args are the parameters, optional */ /* note that since the parameters are input and output, we have to do more work than usual */ for (j=3;j<(MIN(narg,7));j++) { iq = ps[j]; switch (sym[iq].class) { case 8: /* only OK if a single case and the element is in a fp array */ if (sym[iq].type != 3) { printf("GFIT: bad parameter type\n"); printf("subscript of an array not of type FLOAT\n"); return -1; } if (outery != 1) { printf("GFIT: bad parameter type\n"); printf("subscript of an array illegal for multiple fits\n"); return -1; } p = (float *)sym[iq].spec.array.ptr; break; case 4: /* if we need a scalar, just redefine */ if (outery == 1) { redef_scalar(iq,3,&zip); p = &sym[iq].spec.scalar.f; break; } /* if an array, it has to match outery or we redefine */ /* if not float, it gets redefined as float, note that a replace is required since we need to output the same symbol */ if (sym[iq].type != 3) { ana_replace(iq, ana_float(1,&iq) ); } h = (struct ahead *) sym[iq].spec.array.ptr; p = (float *) ((char *)h + sizeof(struct ahead)); nq = h->ndim; n = 1; for(i=0;idims[i]; } if ( outery == n ) { inc = 1; update = 0; } else { if (redef_array(iq, 3, nd-1, &dim[1]) != 1) return -1; h = (struct ahead *) sym[iq].spec.array.ptr; p = (float *) ((char *)h + sizeof(struct ahead)); inc = 1; update = 1; } break; case 1: /* a scalar, if outery is 1, just use pointer, otherwise define an array but save scalar value for first element */ if (sym[iq].type != 3) { ana_replace(iq, ana_float(1,&iq) ); } if (outery == 1) { p = &sym[iq].spec.scalar.f; } else { xq = float_arg(iq); if (redef_array(iq, 3, nd-1, &dim[1]) != 1) return -1; h = (struct ahead *) sym[iq].spec.array.ptr; p = (float *) ((char *)h + sizeof(struct ahead)); *p = xq; inc = 1; update = 1; } break; default: /* define an array or scalar, set to zero */ if (outery == 1) { redef_scalar(iq,3,&zip); p = &sym[iq].spec.scalar.f; } else { if (redef_array(iq, 3, nd-1, &dim[1]) != 1) return -1; ana_zero( 1, &iq); inc = 1; update = 1; } } switch (j) { /* set the real pointer to the temp used above */ case 3: pa = p; inc_a = inc; update_a = update; break; case 4: pb = p; inc_b = inc; update_b = update; break; case 5: px0 = p; inc_x0 = inc; update_x0 = update; break; case 6: pdel = p; inc_del = inc; update_del = update; break; } } /* update flag = 1 only makes sense if outery was > 1, in this case we save the first guess in the guess array. Using a separate guess array allows junmping over bad comvergences when using updates */ if (update_a) guess[0] = *pa; if (update_b) guess[1] = *pb; if (update_x0) guess[2] = *px0; if (update_del) guess[3] = *pdel; tol = 0.005; if (narg > 7) { if (float_arg_stat(ps[7], &tol) != 1) return -1; } niter = 20; /* default iteration count */ if (narg > 8) { if (int_arg_stat(ps[8], &niter) != 1) return -1; } /* for showing progress */ if (narg > 9) { if (int_arg_stat(ps[9], &gfit_show) != 1) return -1; } gfit_converged = 0; /* converged count */ /* ready to call gfit */ while (1) { status = gfit(x,y,f,nxq,pa,pb,px0,pdel,tol,niter); if (status == 1) { gfit_converged += 1; /* update guess only if update flag set and we converged */ if (update_a) guess[0] = *pa; if (update_b) guess[1] = *pb; if (update_x0) guess[2] = *px0; if (update_del) guess[3] = *pdel; } else { if (gfit_warnings) printf("WARNING: GFIT did not converge\n"); } if (--outery <= 0) break; y += nyq; f += nyq; if (outerx > 1) x += nxq; pa += inc_a; pb += inc_b; px0 += inc_x0; pdel += inc_del; if (update_a) *pa = guess[0]; if (update_b) *pb = guess[1]; if (update_x0) *px0 = guess[2]; if (update_del) *pdel = guess[3]; } return 1; } /*------------------------------------------------------------------------- */ int ana_cosinefit(narg,ps) /* the call is COSINEFIT,X,Y,F,[A],[B],[X0],[PERIOD],[TOL],[NITER],[show] */ /* accepts 3 or more arguments */ /* copied from gfit, pretty much the same at this level */ int narg,ps[]; { extern int cosinefit_fix_a, cosinefit_fix_b, cosinefit_fix_x0, cosinefit_fix_del, cosinefit_show; extern int cosinefit_guess_mode, cosinefit_converged, cosinefit_warnings; float *x, *y, *f, guess[4]={0,0,0,0}, tol, zip={0}; /* guess is a,b,x0,del */ float *p, *pa, *pb, *px0, *pdel, xq; int niter, nxq, nd, outerx, j, nyq, outery, dim[8], iq, n, nq, i; int status, inc_a, inc_b, inc_x0, inc_del, inc; int update_a, update_b, update_x0 ,update_del, update; struct ahead *h; iq = ps[0]; if ( sym[iq].class != 4 ) return execute_error(66); /*float x and get pointer */ iq = ana_float(1,&iq); h = (struct ahead *) sym[iq].spec.array.ptr; x = (float *) ((char *)h + sizeof(struct ahead)); nd = h->ndim; nxq = h->dims[0]; outerx = 1; if (nd > 1) for(j=1;jdims[j]; /* now y */ iq = ps[1]; if ( sym[iq].class != 4 ) return execute_error(66); iq = ana_float(1,&iq); h = (struct ahead *) sym[iq].spec.array.ptr; y = (float *) ((char *)h + sizeof(struct ahead)); nd = h->ndim; nyq = h->dims[0]; outery = 1; dim[0] = nyq; if (nd > 1) for(j=1;jdims[j]; dim[j] = h->dims[j]; } /* nxq must = nyq and outerx must = outery or 1 */ if ( nxq != nyq || ( outerx != outery && outerx != 1) ) return execute_error(103); /* now f, make a match for y array */ iq = ps[2]; if (redef_array(iq, 3, nd, dim) != 1) return -1; h = (struct ahead *) sym[iq].spec.array.ptr; f = (float *) ((char *)h + sizeof(struct ahead)); /* check for presence of fixed parameters */ if ( (cosinefit_fix_a && narg < 4) || (cosinefit_fix_b && narg < 5) || (cosinefit_fix_x0 && narg < 6) || (cosinefit_fix_del && narg < 7) ) { printf("cosinefit: missing fixed parameter\n"); return -1; } /* a place to point to for optional args */ pa = &guess[0]; pb = &guess[1]; px0 = &guess[2]; pdel = &guess[3]; inc_a = inc_b = inc_x0 = inc_del = 0; update_a = update_b = update_x0 = update_del = 0; /* the next 4 args are the parameters, optional */ /* note that since the parameters are input and output, we have to do more work than usual */ for (j=3;j<(MIN(narg,7));j++) { iq = ps[j]; switch (sym[iq].class) { case 8: /* only OK if a single case and the element is in a fp array */ if (sym[iq].type != 3) { printf("cosinefit: bad parameter type\n"); printf("subscript of an array not of type FLOAT\n"); return -1; } if (outery != 1) { printf("cosinefit: bad parameter type\n"); printf("subscript of an array illegal for multiple fits\n"); return -1; } p = (float *)sym[iq].spec.array.ptr; break; case 4: /* if we need a scalar, just redefine */ if (outery == 1) { redef_scalar(iq,3,&zip); p = &sym[iq].spec.scalar.f; break; } /* if an array, it has to match outery or we redefine */ /* if not float, it gets redefined as float, note that a replace is required since we need to output the same symbol */ if (sym[iq].type != 3) { ana_replace(iq, ana_float(1,&iq) ); } h = (struct ahead *) sym[iq].spec.array.ptr; p = (float *) ((char *)h + sizeof(struct ahead)); nq = h->ndim; n = 1; for(i=0;idims[i]; } if ( outery == n ) { inc = 1; update = 0; } else { if (redef_array(iq, 3, nd-1, &dim[1]) != 1) return -1; h = (struct ahead *) sym[iq].spec.array.ptr; p = (float *) ((char *)h + sizeof(struct ahead)); inc = 1; update = 1; } break; case 1: /* a scalar, if outery is 1, just use pointer, otherwise define an array but save scalar value for first element */ if (sym[iq].type != 3) { ana_replace(iq, ana_float(1,&iq) ); } if (outery == 1) { p = &sym[iq].spec.scalar.f; } else { xq = float_arg(iq); if (redef_array(iq, 3, nd-1, &dim[1]) != 1) return -1; h = (struct ahead *) sym[iq].spec.array.ptr; p = (float *) ((char *)h + sizeof(struct ahead)); *p = xq; inc = 1; update = 1; } break; default: /* define an array or scalar, set to zero */ if (outery == 1) { redef_scalar(iq,3,&zip); p = &sym[iq].spec.scalar.f; } else { if (redef_array(iq, 3, nd-1, &dim[1]) != 1) return -1; ana_zero( 1, &iq); inc = 1; update = 1; } } switch (j) { /* set the real pointer to the temp used above */ case 3: pa = p; inc_a = inc; update_a = update; break; case 4: pb = p; inc_b = inc; update_b = update; break; case 5: px0 = p; inc_x0 = inc; update_x0 = update; break; case 6: pdel = p; inc_del = inc; update_del = update; break; } } /* update flag = 1 only makes sense if outery was > 1, in this case we save the first guess in the guess array. Using a separate guess array allows junmping over bad comvergences when using updates */ if (update_a) guess[0] = *pa; if (update_b) guess[1] = *pb; if (update_x0) guess[2] = *px0; if (update_del) guess[3] = *pdel; tol = 0.005; if (narg > 7) { if (float_arg_stat(ps[7], &tol) != 1) return -1; } niter = 20; /* default iteration count */ if (narg > 8) { if (int_arg_stat(ps[8], &niter) != 1) return -1; } /* for showing progress */ if (narg > 9) { if (int_arg_stat(ps[9], &cosinefit_show) != 1) return -1; } cosinefit_converged = 0; /* converged count */ /* ready to call cosinefit */ while (1) { status = cosinefit(x,y,f,nxq,pa,pb,px0,pdel,tol,niter); if (status == 1) { cosinefit_converged += 1; /* update guess only if update flag set and we converged */ if (update_a) guess[0] = *pa; if (update_b) guess[1] = *pb; if (update_x0) guess[2] = *px0; if (update_del) guess[3] = *pdel; } else { if (cosinefit_warnings) printf("WARNING: cosinefit did not converge\n"); } if (--outery <= 0) break; y += nyq; f += nyq; if (outerx > 1) x += nxq; pa += inc_a; pb += inc_b; px0 += inc_x0; pdel += inc_del; if (update_a) *pa = guess[0]; if (update_b) *pb = guess[1]; if (update_x0) *px0 = guess[2]; if (update_del) *pdel = guess[3]; } return 1; } /*------------------------------------------------------------------------- */ int ana_lower(narg,ps) /* lower function */ /* convert all letters in a string to lower case */ /* 1/4/98 - accommodate string arrays also */ int narg,ps[]; { byte *p1, *p2; int mq; int result_sym; switch ( sym[ ps[0] ].class ) { case STRING: mq = sym[ps[0]].spec.array.bstore - 1; result_sym = string_scratch(mq); /*for resultant string */ p1 = (byte *) sym[ps[0] ].spec.array.ptr; p2 = (byte *) sym[result_sym].spec.array.ptr; while (mq--) *p2++ = isupper( (int) *p1 ) ? (byte) tolower( (int) *p1++) : *p1++; *p2 = 0; /*ending null */ break; case STRING_PTR: { struct ahead *h; int nd, dim[8], j, n; byte **p, **q; h = (struct ahead *) sym[ps[0]].spec.array.ptr; nd = h->ndim; n = 1; for (j=0;jdims[j]; n *= h->dims[j]; } //printf("nd, n = %d, %d\n", nd, n); q = (byte **) ((byte *) h + sizeof(struct ahead)); result_sym = strarr_scratch(nd, dim); /* we have to loop through all the strings, make a result for each */ h = (struct ahead *) sym[result_sym].spec.array.ptr; p = (byte **) ((byte *) h + sizeof(struct ahead)); for (j=0;jndim; n = 1; for (j=0;jdims[j]; n *= h->dims[j]; } q = (char **) ((char *) h + sizeof(struct ahead)); result_sym = strarr_scratch(nd, dim); /* we have to loop through all the strings, make a result for each */ h = (struct ahead *) sym[result_sym].spec.array.ptr; p = (char **) ((char *) h + sizeof(struct ahead)); for (j=0;jndim; n = 1; for (j=0;jdims[j]; n *= h->dims[j]; } q = (char **) ((char *) h + sizeof(struct ahead)); result_sym = strarr_scratch(nd, dim); /* we have to loop through all the strings, make a result for each */ h = (struct ahead *) sym[result_sym].spec.array.ptr; p = (char **) ((char *) h + sizeof(struct ahead)); for (j=0;j= 4 ) if (int_arg_stat(ps[3], &nc) != 1) return -1; /* do in 2 passes, first count the number of target substrings up to nc */ n= 0; pq = p1; match_size = strlen(p2); rep_size = strlen(p3); while ( (p4 = strstr( p1, p2)) != NULL) { n++; p1 = p4 +match_size; if ( nc > 0 && n >= nc ) break; } /* n is the count to replace, compute storage */ mq = sym[ps[0]].spec.array.bstore - 1 + n * (rep_size - match_size); result_sym = string_scratch(mq); /*for resultant string */ p4 = (char *) sym[ result_sym ].spec.array.ptr; p1 = pq; while ( n-- ) { pq = strstr( p1, p2); /*transfer any before the target */ nc = ( pq - p1 ); while (nc--) *p4++ = *p1++; pq = p3; nc = rep_size; while (nc--) *p4++ = *pq++; p1 += match_size; } /*leftovers */ nc = mq - (p4 - (char *) sym[ result_sym ].spec.array.ptr); if (nc > 0) { while (nc--) { *p4++ = *p1++; } } *p4 = '\0'; return result_sym; } /*------------------------------------------------------------------------- */ int ana_strloc(narg,ps) /* STRLOC function */ /* returns rest of source starting at object string */ int narg,ps[]; { char *p1, *p2, *p3, *p4; int mq, off; int result_sym; string_func_setup /* does p1 and p2 */ p3 = strstr( p1, p2); mq = 0; off = (p3 - p1); if ( p3 != NULL ) mq = sym[ps[0]].spec.array.bstore - 1 - off; result_sym = string_scratch(mq); /*for resultant string */ p4 = (char *) sym[result_sym].spec.array.ptr; if (mq != 0) bcopy( p3, p4, mq ); *(p4 + mq) = 0; return result_sym; } /*------------------------------------------------------------------------- */ int ana_strskp(narg,ps) /* STRSKP function */ /* returns rest of source starting AFTER object string */ int narg,ps[]; { char *p1, *p2, *p3; int mq, off; int result_sym; string_func_setup /* does p1 and p2 */ p3 = strstr( p1, p2); mq = 0; off = (p3 - p1) + sym[ps[1]].spec.array.bstore - 1; if ( p3 != NULL ) mq = sym[ps[0]].spec.array.bstore - 1 - off; if (mq < 0) mq = 0; result_sym = string_scratch(mq); /*for resultant string */ p3 = (char *) sym[result_sym].spec.array.ptr; if (mq != 0) bcopy( p1 + off, p3, mq ); *(p3 + mq) = 0; return result_sym; } /*------------------------------------------------------------------------- */ int ana_strtok(narg,ps) /* STRTOK function */ /* starts a strtok set, basically an interface for the C strtok function the first call needs the input string, after that, specify a NULL string for the first argument to get the next token in that string. The second argument contains the delimiters */ int narg,ps[]; { char *p1, *p2, *p3, *p4; int mq; int result_sym; static char * poriginal; static int size, original_string_sym; if ( sym[ ps[0] ].class != 2 ) return execute_error(70); if ( sym[ ps[1] ].class != 2 ) return execute_error(70); p1 = (char *) sym[ps[0] ].spec.array.ptr; mq = sym[ps[0]].spec.array.bstore; if (mq <= 1) { /* a NULL first arg., we must check that the old one still is reasonable */ p1 = NULL; p2 = (char *) sym[original_string_sym].spec.array.ptr; mq = sym[original_string_sym].spec.array.bstore; if (p2 != poriginal || size != mq) { printf("STRTOK: original string found invalid\n"); return -1; } } else { /* a new start */ poriginal = p1; size = mq; original_string_sym = ps[0]; } p2 = (char *) sym[ps[1] ].spec.array.ptr; p3 = strtok( p1, p2); if ( p3 != NULL ) mq = strlen(p3); else mq = 0; if (mq < 0) mq = 0; result_sym = string_scratch(mq); /*for resultant string */ p4 = (char *) sym[result_sym].spec.array.ptr; if (mq != 0) bcopy( p3, p4, mq ); *(p4 + mq) = 0; return result_sym; } /*------------------------------------------------------------------------- */ int ana_strpbrk(narg,ps) /* SRTPBRK function */ /* returns a string from s1 starting with the first char from s2 found in s1, if none found, an empty string */ int narg,ps[]; { char *p1, *p2, *p3, *p4; int mq, off; int result_sym; string_func_setup /* does p1 and p2 */ p3 = strpbrk( p1, p2); mq = 0; off = (p3 - p1); if ( p3 != NULL ) mq = sym[ps[0]].spec.array.bstore - 1 - off; result_sym = string_scratch(mq); /*for resultant string */ p4 = (char *) sym[result_sym].spec.array.ptr; //printf("off, mq in strpbrk = %d %d\n", off,mq); if (mq != 0) bcopy( p3, p4, mq ); *(p4 + mq) = 0; return result_sym; } /*------------------------------------------------------------------------- */ int ana_skipc(narg,ps) /* SKIPC function */ /* skip over characters */ int narg,ps[]; { char *p1, *p2, *p3; int mq, off; int result_sym; string_func_setup /* does p1 and p2 */ off = strspn( p1, p2); mq = sym[ps[0]].spec.array.bstore - 1 - off; if (mq < 0) mq = 0; result_sym = string_scratch(mq); /*for resultant string */ p3 = (char *) sym[result_sym].spec.array.ptr; if (mq != 0) bcopy( p1 + off, p3, mq ); *(p3 + mq) = 0; return result_sym; } /*------------------------------------------------------------------------- */ int ana_strsub(narg,ps) /* STRSUB function */ /* return a substring */ int narg,ps[]; { char *p1, *p2, *p3; int mq, off, len; int result_sym; if ( sym[ ps[0] ].class != 2 ) return execute_error(70); p1 = (char *) sym[ps[0] ].spec.array.ptr; if (int_arg_stat(ps[1], &off) != 1) return -1; if (int_arg_stat(ps[2], &len) != 1) return -1; if (len < 0 ) len = 0; mq = sym[ps[0]].spec.array.bstore - 1 - off; mq = len <= mq ? len : mq; result_sym = string_scratch(mq); /*for resultant string */ p3 = (char *) sym[result_sym].spec.array.ptr; if (mq != 0) bcopy( p1 + off, p3, mq ); *(p3 + mq) = 0; return result_sym; } /*------------------------------------------------------------------------- */ int ana_strtrim(narg,ps) /* STRTRIM function */ /* trim trailing blanks and tabs and any control chars */ int narg,ps[]; { /* 9/27/97 - modified to return a zero length if all are blanks, etc */ char *p1, *p2, *p3; int mq, n; int result_sym; if ( sym[ ps[0] ].class != 2 ) return execute_error(70); p1 = (char *) sym[ps[0] ].spec.array.ptr; mq = sym[ps[0]].spec.array.bstore - 1; /* length of input string */ if (mq != 0) { p2 = p1 + mq; /*just after end */ /* while (1) { --p2; if ( *p2 > 32 || p2 <= p1 ) break; } */ /* 9/27/97 mod */ while (1) { --p2; if ( p2 < p1 || *p2 > 32 ) break; } mq = p2 - p1 + 1; } result_sym = string_scratch(mq); /*for resultant string */ p3 = (char *) sym[result_sym].spec.array.ptr; if (mq != 0) bcopy( p1, p3, mq ); *(p3 + mq) = 0; return result_sym; } /*------------------------------------------------------------------------- */ int ana_istring(narg,ps) /* ISTRING function */ /* string in I format, 3 arguments with last 2 optional, call is s = istring(value, width, mode) where width is the field size (defaults is 10), and mode=0 => leading blanks, 1 => leading 0's, 2 => width = # digits */ int narg,ps[]; { char s[6], *p3; int iq, mq, n, k, mode; int result_sym; if (int_arg_stat(ps[0], &k) != 1) return -1; n=10; mode = 0; if ( narg > 1 ) if (int_arg_stat(ps[1], &n) != 1) return -1; if ( narg > 2 ) if (int_arg_stat(ps[2], &mode) != 1) return -1; switch (mode) { default: case 0: sprintf(s, "%%%1dd", n); break; case 1: sprintf(s, "%%0%1dd", n); break; case 2: sprintf(s, "%%1d"); break; } /* s now has the control string for the formatting, use line to store result since length can be uncertain for some cases */ sprintf(line, s, k); mq = strlen(line); /*printf("mq = %d, s = %s, line = %s\n",mq,s,line);*/ result_sym = string_scratch(mq); /*for resultant string */ p3 = (char *) sym[result_sym].spec.array.ptr; bcopy( line, p3, mq); *(p3 + mq) = 0; return result_sym; } /*------------------------------------------------------------------------- */ int ana_hexstring(narg,ps) /* HEXSTRING function */ /* an integer in hex format */ int narg,ps[]; { char s[6], *p3; int iq, mq, n, k, mode; int result_sym; if (int_arg_stat(ps[0], &k) != 1) return -1; n=1; mode = 0; if ( narg > 1 ) if (int_arg_stat(ps[1], &n) != 1) return -1; if ( narg > 2 ) if (int_arg_stat(ps[2], &mode) != 1) return -1; switch (mode) { default: case 0: sprintf(s, "%%#%1dx", n); break; case 1: sprintf(s, "%%%1dx", n); break; } /* s now has the control string for the formatting, use line to store result since length can be uncertain for some cases */ sprintf(line, s, k); mq = strlen(line); /*printf("mq = %d, s = %s, line = %s\n",mq,s,line);*/ result_sym = string_scratch(mq); /*for resultant string */ p3 = (char *) sym[result_sym].spec.array.ptr; bcopy( line, p3, mq); *(p3 + mq) = 0; return result_sym; } /*------------------------------------------------------------------------- */ int ana_fstring(narg,ps) /* FSTRING function */ /* formatted fp string, limited */ int narg,ps[]; { char s[6], *p3, *p1; int iq, mq, n, mode; int result_sym; float xf; if ( sym[ ps[0] ].class != 2 ) return execute_error(70); p1 = (char *) sym[ps[0]].spec.array.ptr; /* p1 is the format string */ iq = ps[1]; if (float_arg_stat(iq, &xf) != 1) return -1; /* p1 now has the control string for the formatting, use line to store result since length can be uncertain for some cases */ sprintf(line, p1, xf); mq = strlen(line); /*printf("mq = %d, s = %s, line = %s\n",mq,s,line);*/ result_sym = string_scratch(mq); /*for resultant string */ p3 = (char *) sym[result_sym].spec.array.ptr; bcopy( line, p3, mq); *(p3 + mq) = 0; return result_sym; } /*------------------------------------------------------------------------- */ int ana_ifstring(narg,ps) /* IFSTRING function */ /* formatted int string, limited */ int narg,ps[]; { char s[6], *p3, *p1; int iq, mq, n, mode; int result_sym; int k; if ( sym[ ps[0] ].class != 2 ) return execute_error(70); p1 = (char *) sym[ps[0]].spec.array.ptr; /* p1 is the format string */ iq = ps[1]; if (int_arg_stat(iq, &k) != 1) return -1; /* p1 now has the control string for the formatting, use line to store result since length can be uncertain for some cases */ sprintf(line, p1, k); mq = strlen(line); result_sym = string_scratch(mq); /*for resultant string */ p3 = (char *) sym[result_sym].spec.array.ptr; bcopy( line, p3, mq); *(p3 + mq) = 0; return result_sym; } /*------------------------------------------------------------------------- */ int ana_sscanf(int narg,int ps[]) /* SSCANF function */ /* interface to C sscanf, limited to 31 values here */ { extern int scrat[NSCRAT]; char *format, *s, *sstart, C, F, *fstart, save; int nc, *conversions, nconv, isym, iq, itype, otype, scalarFlag, in1, in2; char formbuf[32]; union scanresult { byte b; short w; int i; float f; double d; long int li; long long int lli; long double ld;} scanresult; union types_ptr p1; union scalar xs, *xsp; void *pscanresult = &scanresult.b; int NoAssign, IsShort, IsLongDoub, IsLong, HaveWidth, Width, nconsumed; if ( sym[ ps[0] ].class != 2 ) return execute_error(70); s = sstart = (char *) sym[ps[0]].spec.array.ptr; /* s in the input string */ if ( sym[ ps[1] ].class != 2 ) return execute_error(70); format = (char *) sym[ps[1]].spec.array.ptr; /* format is the format string */ result_sym = scalar_scratch(2); /* for the return, this will be # of conversions */ conversions = &sym[result_sym].spec.scalar.l; *conversions = 0; nconsumed = 0; /* get the first output argument, has to be at least one */ isym = 2; /* this is loosely based on the FreeBSD _scanf, we have to parse the formats in order to set up the pointers properly for the ana variables which include strings, but we use the system sscanf for anything hard in a single conversion call */ /* we have to avoid collisions between input and output variables */ in1 = ps[0]; in2 = ps[1]; /* Walk over the format string */ while (F = *format++) { /* anything left in the input string ? */ if (*s == '\0') goto WrapUP; /* Check for a conversion */ if (F != '%' || *format == '%') { /* %% or any char other than % */ if (F == '%') { ++format; } /* Check for a match */ if (isspace (F)) { /* Special white space handling: Any whitespace matches * any amount of whitespace including none(!). So this * match will never fail. */ while (isspace (*s)) { s++; } continue; } else if (F != *s) { /* A mismatch. We will stop scanning the input and return * the number of conversions. */ goto WrapUP; } else { /* A match. Read the next input character and start over */ s++; continue; } } else { /* A conversion. Skip the percent sign. */ fstart = format - 1; /* keep this positon in the format string*/ F = *format++; /* we are going to let the real sscanf do each individual variable, but we have to span over the format specifier and setup the appropiate ana variable */ /* Initialize variables */ NoAssign = 0; IsShort = 0; IsLongDoub = 0; IsLong = 0; HaveWidth = 0; Width = 4*NSCRAT; /* Check for flags. */ while (1) { if (isdigit (F)) { HaveWidth = 1; Width = 0; do { /* don't need the width here except for string format */ /* ### Non portable ### */ Width = Width * 10 + (F & 0x0F); F = *format++; } while (isdigit (F)); } else { switch (F) { case '*': NoAssign = 1; break; case 'h': IsShort += 1; if (IsShort > 2) IsShort = 2; break; case 'l': IsLong += 1; if (IsLong > 2) IsLong = 2; break; case 'L': IsLongDoub = 1; break; default: goto FlagsDone; } F = *format++; } } FlagsDone: scalarFlag = 0; /* Check for the actual conversion character */ switch (F) { case 'D': printf("The %%D format is obsolete, use %%ld instead"); return result_sym; break; case 'u': case 'd': case 'i': case 'o': case 'x': case 'X': case 'E': case 'e': case 'f': case 'g': nc = format - fstart; //printf("nc = %d\n", nc); save = *format; *format = 0; snprintf(formbuf, 30, "%*s%%n", nc, fstart); //printf("formbuf is: %s\n", formbuf); /* all these will put a scalar value in pscanresult */ nconv = sscanf(s, formbuf, pscanresult, &nc); *format = save; *conversions += 1; //printf("nc = %d, conversions = %d\n", nc, *conversions); s = s + nc; /* our scalar is in pscanresult, we sort it out in another switch statement below */ scalarFlag = 1; break; case 's': /* Whitespace terminated string */ C = *s; nc = 0; while (isspace (C)) { C = *s++; } if (!NoAssign) { char *p, *sout; int n; sout = (char *) scrat; while (!isspace (C) && Width--) { *sout++ = C; s++; C = *s; nc++; if (C == '\0') break; } /* copy the string just read to an ana string */ iq = ps[isym++]; if (iq == in1 || iq == in2) return execute_error(138); redef_string (iq, nc ); p = (char *) sym[iq].spec.array.ptr; sout = (char*) scrat; n = nc; while (n--) *p++ = *sout++; *p = 0; /*load it */ } else { /* just scan until a null or EOF, but still a limit */ while (!isspace (C) && Width--) { C = *s++; } } *conversions += 1; break; case 'c': /* Fixed length string, NOT zero terminated */ C = *s; nc = 0; if (!HaveWidth) { /* No width given, default is 1 */ Width = 1; } //printf("c, width = %d, s = %#x, value in C = %d\n", Width, s, C); if (!NoAssign) { char *p, *sout; int n; sout = (char *) scrat; while (Width--) { *sout++ = C; s++; C = *s; nc++; //printf("new value in C = %d\n", C); if (C == '\0') break; } //printf("c, nc = %d\n", nc); //printf("s = %#x, string: %s\n", s, s); /* copy the string just read to an ana string */ iq = ps[isym++]; if (iq == in1 || iq == in2) return execute_error(138); redef_string (iq, nc ); p = (char *) sym[iq].spec.array.ptr; sout = (char*) scrat; n = nc; while (n--) *p++ = *sout++; *p = 0; /*load it */ } else { /* Just skip as many chars as given */ while (Width--) { C = *s++; } } *conversions += 1; //printf("end s = %#x, string: %s\n", s, s); break; /* the [ format looks like more work and hasn't been done yet */ // case '[': // /* String using characters from a set */ // Invert = 0; // /* Clear the set */ // memset (CharSet, 0, sizeof (CharSet)); // F = *format++; // if (F == '^') { // Invert = 1; // F = *format++; // } // if (F == ']') { // AddCharToSet (']'); // F = *format++; // } // /* Read the characters that are part of the set */ // while (F != ']' && F != '\0') { // if (*format == '-') { // /* A range. Get start and end, skip the '-' */ // Start = F; // F = *++format; // ++format; // if (F == ']') { // /* '-' as last char means: include '-' */ // AddCharToSet (Start); // AddCharToSet ('-'); // } else if (F != '\0') { // /* Include all chars in the range */ // while (1) { // AddCharToSet (Start); // if (Start == F) { // break; // } // ++Start; // } // /* Get next char after range */ // F = *format++; // } // } else { // /* Just a character */ // AddCharToSet (F); // /* Get next char */ // F = *format++; // } // } // // /* Invert the set if requested */ // if (Invert) { // for (Start = 0; Start < sizeof (CharSet); ++Start) { // CharSet[Start] ^= 0xFF; // } // } // // /* We have the set in CharSet. Read characters and // * store them into a string while they are part of // * the set. // */ // if (!NoAssign) { // S = va_arg (ap, char*); // while (IsCharInSet (C) && Width--) { // *S++ = C; // ReadChar (); // } // *S = '\0'; // } else { // while (IsCharInSet (C) && Width--) { // ReadChar (); // } // } // ++Conversions; // break; // /* note that %p doesn't make sense for ana */ case 'n': /* Store characters consumed so far */ nc = s - sstart; *conversions += 1; scalarFlag = 1; break; default: /* Invalid conversion */ printf("sscanf - Invalid or unsupported conversion %c\n", F); return result_sym; break; } if (scalarFlag) { if (NoAssign == 0) { /* and check our format again, scalars are sent to ana variables */ switch (F) { /* case D would have caused an exit above */ /* these are int's unless we had a flag set for short, char, long, or long long, we put everything into an int (ana doesn't do 64 or 132 bit integers) but don't extend the sign for the unsigned cases uox if they have a h or hh, the l and ll cases suffer truncation */ case 'u': case 'o': case 'x': case 'X': /* these are unsigned, */ if (IsShort) { if (IsShort == 2) xs.l = ((int) scanresult.b) & 0x0f; else xs.l = ((int) scanresult.w) & 0xff; } else if (IsLong) { if (IsLong == 2) xs.l = (int) scanresult.lli; else xs.l = (int) scanresult.li; } else { xs.l = scanresult.i; } //printf("unsigned result = %d\n", xs.l); itype = 2; break; case 'd': case 'i': /* these are signed so we sign extend the short ones, the int cast should do this */ if (IsShort) { if (IsShort == 2) xs.l = (int) scanresult.b; else xs.l = (int) scanresult.w; } else if (IsLong) { if (IsLong == 2) xs.l = (int) scanresult.lli; else xs.l = (int) scanresult.li; } else { xs.l = scanresult.i; } //printf("signed result = %d\n", xs.l); itype = 2; break; case 'n': xs.l = nc; itype = 2; break; case 'E': case 'e': case 'f': case 'g': /* these are all the same but there might have been an l or L modifier for any of them */ if (IsLong) { xs.d = scanresult.d; } else if (IsLongDoub) { xs.d = (double) scanresult.ld; /* generally truncated */ } else { xs.d = (double) scanresult.f; } //printf("double result = %g\n", xs.d); itype = 4; break; } if (isym == in1 || isym == in2) return execute_error(138); iq = ps[isym++]; xsp = &xs; //printf("iq = %d, sym[iq].class = %d\n", isym, sym[iq].class); /* only scalars are supported for non-string types */ switch ( sym[iq].class ) { case 1: p1.l = &sym[iq].spec.scalar.l; ana_array_convert( &xsp, &p1.l, itype, sym[iq].type, 1); break; case 8: p1.l = sym[iq].spec.array.ptr; ana_array_convert( &xsp, &p1.l, itype, sym[iq].type, 1); break; case 255: /* undefined, so match itype */ if (redef_scalar(iq, itype, xsp) != 1) return -1; break; default: /* anything else (like an array) is not legal here, strings are done in a different section */ return execute_error(137); } } } if (*s == '\0') goto WrapUP; } } WrapUP: /* Return the number of conversions */ return result_sym; } /*------------------------------------------------------------------------- */ int ana_sprintf(narg,ps) /* SPRINTF function */ /* interface to C sprintf, limited to 31 values here */ int narg,ps[]; { char *p3, *p, *fmt, *fmt2, *p2, *out, save; int iq, mq, n, mode, total_width, fmtend, ns, i, star_flag = 0; int result_sym; int k, k1, k2; double zq; if ( sym[ ps[0] ].class != 2 ) return execute_error(70); p = fmt = (char *) sym[ps[0]].spec.array.ptr; /* p is the format string */ /* get the probable length of the result, may need some parameters */ n = 0; /* next arg, skipping the format string */ /* cribbed from the vasprintf code in the lib library */ /* the estimate is usually quite excessive but safe */ total_width = fmtend = strlen (fmt) + 1; while (*p != '\0') { if (*p++ == '%') { n++; while (strchr ("-+ #0", *p)) ++p; if (*p == '*') { /* an annoying *, we need the parameter now */ ++p; iq = ps[n]; if (n++ >= narg) return execute_error(130); if (int_arg_stat(iq, &k) != 1) return execute_error(129); total_width += abs(k); } else total_width += strtoul (p, (char **)&p, 10); if (*p == '.') { ++p; if (*p == '*') { ++p; /* another *, need the corresponding variable */ iq = ps[n]; if (n++ >= narg) return execute_error(130); if (int_arg_stat(iq, &k) != 1) return execute_error(129); total_width += abs(k); } else total_width += strtoul (p, (char **)&p, 10); } while (strchr ("hlL", *p)) ++p; /* and add some for the defaults and possible expansions */ total_width += 30; switch (*p) { case 'd': case 'i': case 'o': case 'u': case 'x': case 'X': case 'c': case 'f': case 'e': case 'E': case 'g': case 'G': case 'p': case 'n': break; case 's': //printf("a string, n = %d\n", n); if ( sym[ ps[n] ].class != 2 ) return execute_error(70); total_width += strlen ((char *) sym[ps[n]].spec.array.ptr); break; default: printf("SPRINTF - illegal format specifer %c\n", *p); return execute_error(129); break; } } } //printf("size = %d\n", total_width); /* now we malloc this much for string, note this already has an extra for the null from the +1 at the very beginning*/ result_sym = string_scratch(total_width); /*for resultant string */ out = (char *) sym[result_sym].spec.array.ptr; /* and now do the loop again for real */ p = fmt; fmt2 = p; i = 0; /* position in output string */ n = 1; /* next arg */ while (*p != '\0') { if (*p++ == '%') { while (strchr ("-+ #0", *p)) ++p; star_flag = 0; if (*p == '*') { /* an annoying *, we need the parameter now */ ++p; iq = ps[n]; if (n++ >= narg) return execute_error(130); if (int_arg_stat(iq, &k1) != 1) return execute_error(129); star_flag++; } else p = p + strspn(p, "0123456789"); if (*p == '.') { ++p; if (*p == '*') { ++p; /* another *, need the corresponding variable */ iq = ps[n]; if (n++ >= narg) return execute_error(130); if (int_arg_stat(iq, &k2) != 1) return execute_error(129); if (!star_flag) k1 = k2; /* if this is the only one, use k1 */ star_flag++; } else p = p + strspn(p, "0123456789"); } while (strchr ("hlL", *p)) ++p; iq = ps[n]; if (n++ >= narg) return execute_error(130); switch (*p) { case 'd': case 'i': case 'o': case 'u': case 'x': case 'X': case 'c': /* these are all integers */ if (int_arg_stat(iq, &k) != 1) return execute_error(131); /* get the start of the format string and terminate it with a temporary null at the appropiate place, save that character */ save = p[1]; p[1] = 0; switch (star_flag) { case 0: ns = sprintf( out+i, fmt2, k); break; case 1: ns = sprintf( out+i, fmt2, k1, k); break; case 2: ns = sprintf( out+i, fmt2, k1, k2, k); break; } break; case 'f': case 'e': case 'E': case 'g': case 'G': /* these will all work with a double fp scalar */ if (double_arg_stat(iq, &zq) != 1) return execute_error(131); save = p[1]; p[1] = 0; switch (star_flag) { case 0: ns = sprintf( out+i, fmt2, zq); break; case 1: ns = sprintf( out+i, fmt2, k1, zq); break; case 2: ns = sprintf( out+i, fmt2, k1, k2, zq); break; } break; case 'p': /* pointer case, print the location of the ana symbol ? */ case 'n': printf("p and n formats not supported yet\n"); return -1; break; case 's': if (sym[iq].class != 2 ) return execute_error(70); p2 = (char *) sym[iq].spec.array.ptr; save = p[1]; p[1] = 0; switch (star_flag) { case 0: ns = sprintf( out+i, fmt2, p2); break; case 1: ns = sprintf( out+i, fmt2, k1, p2); break; case 2: ns = sprintf( out+i, fmt2, k1, k2, p2); break; } break; default: printf("SPRINTF - illegal format specifer %c\n", *p); /* but restore the format string before we return */ if (*p) p[1] = save; return -1; break; } p[1] = save; fmt2 = ++p; i = i + ns; } } /* anything left ?, if fmt2 is not equal to p, do the rest */ if (fmt2 < p) { ns = sprintf( out+i, fmt2); i = i + ns; } /* terminate the ana string at the amount actually used, don't copy though this means there will be some extra malloced which free will take care of when the string is deleted */ sym[result_sym].spec.array.bstore = i + 1; //printf("i, strlen(out) = %d %d\n", i, strlen(out)); return result_sym; } /*------------------------------------------------------------------------- */ int ana_string(narg,ps) /* STRING function */ /* string from scalar, no format control yet */ int narg,ps[]; { char s[6], *p3; int iq, mq, n, k, mode; int result_sym; float x; double xx; iq = ps[0]; switch (sym[iq].class) { default: return execute_error(10); case 8: iq = class8_to_1(iq); /*scalar ptr case */ case 1: /*scalar case */ switch (sym[iq].type) { case 0: case 1: case 2: k = int_arg(iq); sprintf(line, "%d", k); break; case 3: x = float_arg(iq); sprintf(line, "%g", x); break; case 4: xx = double_arg(iq); sprintf(line, "%g", xx); break; } } mq = strlen(line); /*printf("mq = %d, s = %s, line = %s\n",mq,s,line);*/ result_sym = string_scratch(mq); /*for resultant string */ p3 = (char *) sym[result_sym].spec.array.ptr; bcopy( line, p3, mq); *(p3 + mq) = 0; return result_sym; } /*------------------------------------------------------------------------- */ int ana_strlen(narg,ps) /* strlen function */ /* length, not counting null; duplicates num_elem but insists on string */ int narg,ps[]; { int result_sym; if ( sym[ ps[0] ].class != 2 ) return execute_error(70); result_sym = scalar_scratch(2); sym[result_sym].spec.scalar.l = sym[ps[0]].spec.array.bstore - 1; return result_sym; } /*------------------------------------------------------------------------- */ int ana_decomp(narg,ps) /* decomp routine */ /*replaces a square matrix with it's lu decomposition */ int narg,ps[]; { int iq, nd, nx; struct ahead *h; union types_ptr q1; iq = ps[0]; if ( iq >= temp_base ) return execute_error(104); if ( sym[iq].class != 4 ) return execute_error(66); h = (struct ahead *) sym[iq].spec.array.ptr; nd = h->ndim; nx = h->dims[0]; if ( nd != 2 || nx != h->dims[1] ) return execute_error(76);/*must be 2-D */ /* and square */ if (sym[iq].type < 4 ) iq = ana_float(1, &iq); if (iq != ps[0] ) ana_replace(ps[0], iq); iq = ps[0]; h = (struct ahead *) sym[iq].spec.array.ptr; q1.f = (float *) ((char *)h + sizeof(struct ahead)); /*could be either float or double, support both */ switch ( sym[iq].type ) { case 3: f_decomp( q1.f, nx, nx ); break; case 4: d_decomp( q1.d, nx, nx ); break; } return 1; } /*------------------------------------------------------------------------- */ int ana_dsolve(narg,ps) /* dsolve routine */ /*solves a linear system given lu decomp and rhs */ int narg,ps[]; { int iq, jq, nd, nx, toptype, j, outer; struct ahead *h; union types_ptr q1, q2; iq = ps[0]; /*the lu decomp matrix */ if ( sym[iq].class != 4 ) return execute_error(66); toptype = sym[iq].type; h = (struct ahead *) sym[iq].spec.array.ptr; nd = h->ndim; nx = h->dims[0]; if ( nd != 2 || nx != h->dims[1] ) return execute_error(76);/*must be 2-D */ /* and square */ jq = ps[1]; /*rhs side (may be more than one)*/ if ( jq >= temp_base ) return execute_error(104); toptype = toptype > sym[jq].type ? toptype : sym[iq].type; if (toptype < 3) toptype = 3; h = (struct ahead *) sym[jq].spec.array.ptr; nd = h->ndim; if ( h->dims[0] != nx) return execute_error(79); outer = 1; if (nd > 1) for(j=1;jdims[j]; /*check type and upgrade as necessary */ switch (toptype) { case 3: iq = ana_float( 1, &iq); jq = ana_float( 1, &jq); break; case 4: iq = ana_double( 1, &iq); jq = ana_double( 1, &jq); break; } if (jq != ps[1] ) ana_replace(ps[1], jq); jq = ps[1]; h = (struct ahead *) sym[iq].spec.array.ptr; q1.l = (int *) ((char *)h + sizeof(struct ahead)); h = (struct ahead *) sym[jq].spec.array.ptr; q2.l = (int *) ((char *)h + sizeof(struct ahead)); while (outer--) { /*could be either float or double, support both */ switch ( toptype ) { case 3: f_solve( q1.f, q2.f, nx, nx ); q2.f += nx; break; case 4: d_solve( q1.d, q2.d, nx, nx ); q2.d += nx; break; } } return 1; } /*------------------------------------------------------------------------- */ int ana_matmul(narg,ps) /* matmul function */ /* multiply 2 matrices, everybody upgraded to F*4 or F*8 */ int narg,ps[]; { int iq, jq, nd, nx, toptype, j, outer, ny, nx2, ny2, dim[2], result_sym; struct ahead *h; union types_ptr q1, q2, q3; iq = ps[0]; if ( sym[iq].class != 4 ) return execute_error(66); toptype = sym[iq].type; h = (struct ahead *) sym[iq].spec.array.ptr; nd = h->ndim; nx = h->dims[0]; if ( nd != 2) return execute_error(101);/*must be 2-D */ ny = h->dims[1]; jq = ps[1]; /* second matrix */ toptype = toptype > sym[jq].type ? toptype : sym[jq].type; if (toptype < 3) toptype = 3; h = (struct ahead *) sym[jq].spec.array.ptr; nd = h->ndim; nx2 = h->dims[0]; if ( nd > 2) return execute_error(92);/*must be 1 or 2-D */ if (nd == 1) ny2 = 1; else ny2 = h->dims[1]; /* the second dimension of array 1 must match first of array 2 */ /* if an inner product of 2 vectors is desired, they must be dimensioned as (1,n) and (n,1) and the result is (1,1) */ if (ny != nx2) return execute_error(77); dim[0] = nx; dim[1] = ny2; /* note that if array 2 is a vector then nd = 1 and the result is also a vector */ /*check type and upgrade as necessary */ switch (toptype) { case 3: result_sym = array_scratch(3, nd, dim); iq = ana_float( 1, &iq); jq = ana_float( 1, &jq); break; case 4: result_sym = array_scratch(4, nd, dim); iq = ana_double( 1, &iq); jq = ana_double( 1, &jq); break; } h = (struct ahead *) sym[iq].spec.array.ptr; q1.l = (int *) ((char *)h + sizeof(struct ahead)); h = (struct ahead *) sym[jq].spec.array.ptr; q2.l = (int *) ((char *)h + sizeof(struct ahead)); h = (struct ahead *) sym[result_sym].spec.array.ptr; q3.l = (int *) ((char *)h + sizeof(struct ahead)); /*could be either float or double, support both */ switch ( toptype ) { register int n1, n2, n3; case 3: { register float sum; register float *a, *b, *c, *abase, *bbase; abase = q3.f; bbase = q1.f; n1 = nx; while (n1--) { n3 = ny2; a = abase++; c = q2.f; while (n3--) { sum = 0.; n2 = ny; b = bbase; while (n2--) { sum += *b * *c++; b += nx; } *a = sum; a += nx; } bbase++; } } break; case 4: { register double sum; register double *a, *b, *c, *abase, *bbase; abase = q3.d; bbase = q1.d; n1 = nx; while (n1--) { n3 = ny2; a = abase++; c = q2.d; while (n3--) { sum = 0.; n2 = ny; b = bbase; while (n2--) {sum += *b * *c++; b += nx; } *a = sum; a += nx; } bbase++; } } break; // #if NeXT // matmul(q3.f, q1.f, q2.f, &nx, &ny,&ny2); break; // case 4: // matmuld(q3.d, q1.d, q2.d, &nx, &ny, &ny2); break; // } // #else // matmul_(q3.f, q1.f, q2.f, &nx, &ny, &ny2); break; // case 4: // { register double sum; // matmuld_(q3.d, q1.d, q2.d, &nx, &ny,&ny2); break; // } // #endif } return result_sym; } /*------------------------------------------------------------------------- */ int ana_pit(narg,ps) /* pit function */ /* polynomial fit, CALL IS C=PIT([X],Y,[NPOW])*/ int narg,ps[]; { double *a, *fbase, *cfbase; union types_ptr qxbase,qybase; int npow, symx, symy, nd, nxq, outer, outerx, iq, jq, dim[2], cfsym; int j, i, n, ib, ie, toptype, k; struct ahead *h; /* depending on # of nargs, we may have to default x or npow, narg=2 case is ambiguous until we get class */ npow = 0; symx = -1; switch (narg) { case 3: symx = ps[0]; symy = ps[1]; npow = int_arg( ps[2] ); break; case 2: /* two possibilities, determined by class of ps[1] */ if ( sym[ps[1]].class == 4 ) { /* the (x,y) case */ symx = ps[0]; symy = ps[1]; } else { symy = ps[0]; npow = int_arg( ps[1] ); } break; case 1: /* just the y vector */ symy = ps[0]; break; } if (npow < 1 || npow >10 ) return execute_error(78); /*check power */ if (sym[symy].class != 4 ) return execute_error(66); /* if either y or x (if x specified) is double, do a double calc. */ toptype = sym[symy].type < 3 ? 3 : sym[symy].type; h = (struct ahead *) sym[symy].spec.array.ptr; /*array may be upgraded but y dimensions will be same, get now */ nd = h->ndim; nxq = h->dims[0]; outer = 1; if (nd > 1) for(j=1;jdims[j]; /* check x if specified */ if (symx > 0) { if (sym[symx].class != 4 ) return execute_error(66); toptype = sym[symx].type > toptype ? sym[symx].type : toptype; h = (struct ahead *) sym[symx].spec.array.ptr; /* a specified x must match y in some regards */ nd = h->ndim; if (h->dims[0] != nxq) return execute_error(77); outerx = 1; if (nd > 1) for(j=1;jdims[j]; if (outerx != 1 && outerx != outer ) return execute_error(77); } else { /* make a fake x array */ symx = setxpit(toptype, nxq); outerx = 1; } /* check types and upgrade as necessary */ switch (toptype) { case 3: symy = ana_float( 1, &symy); symx = ana_float( 1, &symx); break; case 4: symy = ana_double( 1, &symy); symx = ana_double( 1, &symx); break; } h = (struct ahead *) sym[symy].spec.array.ptr; qybase.l = (int *) ((char *)h + sizeof(struct ahead)); h = (struct ahead *) sym[symx].spec.array.ptr; qxbase.l = (int *) ((char *)h + sizeof(struct ahead)); /* get the various matrices and vectors needed */ dim[0] = npow + 1; dim[1] = outer; if (outer > 1) cfsym = array_scratch(4, 2, dim); /*always double */ else cfsym = array_scratch(4, 1, dim); h = (struct ahead *) sym[cfsym].spec.array.ptr; cfbase = (double*) ((char *)h + sizeof(struct ahead)); /*scratch array */ if ( (fbase = (double *) malloc( nxq * sizeof(double))) == NULL) return execute_error(23); if ( (a = (double *) malloc((npow+1)*(npow+1) * sizeof(double))) == NULL) return execute_error(23); /*loop over the number of individual fits to do */ while (outer--) { /* ready to start actual calculation, all done in cana_pit */ cana_pit( &qxbase, &qybase, nxq, npow, a, cfbase, fbase, toptype, outerx); cfbase += npow+1; } /*end of while (outer--) loop */ free(a); free(fbase); return cfsym; } /*------------------------------------------------------------------------- */ int cana_pit( qxbase, qybase, nxq, npow, a, cfbase, fbase, toptype, outerx ) /* internal routine used by ana_pit, ana_trend, and ana_detrend */ union types_ptr *qxbase,*qybase; double *a, *fbase, *cfbase; int nxq, npow, toptype, outerx; { double sum, *f, *cf; union types_ptr qx,qy; int i, n, ib, ie, k; cf = cfbase; f = fbase; qy.l = qybase->l; qx.l = qxbase->l; switch (toptype) { case 3: /* printf("float case\n"); for (k=0;kl; while (n--) {sum += *f; *f = (*f) * (*qx.f++); f++; } ib = (k+1-npow) > 0 ? (k+1-npow) : 0; /*max of */ ie = (k+1) > (npow) ? (npow) : (k+1); /*min of */ for (i=ib;i<=ie;i++) { *( a + i + (npow+1)*(ie+ib-i) ) = sum; } } /* now the rhs */ f = fbase; n = nxq; while (n--) *f++ = *qy.f++; for (k=0;k <= (npow);k++) { sum = 0.0; f = fbase; n = nxq; qx.l = qxbase->l; while (n--) {sum += *f; *f = (*f) * (*qx.f++); f++; } *cf++ = sum; } qybase->f += nxq; if (outerx != 1) qxbase->f += nxq; break; case 4: /* printf("double case\n"); for (k=0;kl; while (n--) {sum += *f; *f = (*f) * (*qx.d++); f++; } ib = (k+1-npow) > 0 ? (k+1-npow) : 0; /*max of */ ie = (k+1) > (npow) ? (npow) : (k+1); /*min of */ for (i=ib;i<=ie;i++) { *( a + i + (npow+1)*(ie+ib-i) ) = sum; } } /* now the rhs */ f = fbase; n = nxq; while (n--) *f++ = *qy.d++; for (k=0;k <= (npow);k++) { sum = 0.0; f = fbase; n = nxq; qx.l = qxbase->l; while (n--) {sum += *f; *f = (*f) * (*qx.d++); f++; } *cf++ = sum; } qybase->d += nxq; if (outerx != 1) qxbase->d += nxq; break; } d_decomp( a, npow+1, npow+1); d_solve(a, cfbase, npow+1, npow+1); /*for (k=0;k<=npow;k++) printf("k, cf(k) %d %f\n",k,*(cfbase+k));*/ return 1; } /*------------------------------------------------------------------------- */ int setxpit(type,n) /*used by several routines to setup an x array */ /*consisting of {0,1/n,2/n,...,(n-1)/n} */ /* type must be 3 or 4 */ int type,n; { union types_ptr qx; float del, tmp; double ddel, dtmp; int nsym, i, nx; struct ahead *h; nx = n; nsym = array_scratch(type, 1, &nx ); h = (struct ahead *) sym[nsym].spec.array.ptr; qx.l = (int *) ((char *)h + sizeof(struct ahead)); switch (type) { case 3: del = 1.0 / nx; *qx.f = 0.0; while (--nx) *(qx.f+1) = *qx.f++ + del; break; case 4: ddel = 1.0/ nx; *qx.d = 0.0; while (--nx) *(qx.d+1) = *qx.d++ + ddel; break; } qx.l = (int *) ((char *)h + sizeof(struct ahead)); return nsym; } /*------------------------------------------------------------------------- */ int ana_detrend(narg,ps) /*detrend function */ /* detrend using a polynomial fit, CALL IS DT=detrend(Y,[NPOW])*/ int narg,ps[]; { /* just call trend with the detrend flag set */ detrend_flag = 1; return ana_trend(narg,ps); } /*------------------------------------------------------------------------- */ int ana_trend(narg,ps) /*trend function */ /* trend using a polynomial fit, CALL IS T=trend(Y,[NPOW])*/ int narg,ps[]; { union types_ptr qx,qy; double sum, *a, *f, *cf; double *abase, *fbase, *cfbase; union types_ptr qxbase,qybase,qzbase; int npow, symx, symy, nd, nxq, outer, outerx, iq, jq, dim[2]; int j, i, n, ib, ie, toptype, k, result_sym; struct ahead *h; /* second argument is optional, default is linear fit (npow = 1) */ npow = 1; switch (narg) { case 2: npow = int_arg( ps[1] ); case 1: symy = ps[0]; break; } if (npow < 1 || npow >10 ) return execute_error(78); /*check power */ if (sym[symy].class != 4 ) return execute_error(66); toptype = sym[symy].type < 3 ? 3 : sym[symy].type; h = (struct ahead *) sym[symy].spec.array.ptr; /*array may be upgraded but y dimensions will be same, get now */ nd = h->ndim; nxq = h->dims[0]; outer = 1; if (nd > 1) for(j=1;jdims[j]; /* make a fake x array */ symx = setxpit(toptype, nxq); outerx = 1; /* check type and upgrade as necessary */ switch (toptype) { case 3: symy = ana_float( 1, &symy); symx = ana_float( 1, &symx); break; case 4: symy = ana_double( 1, &symy); symx = ana_double( 1, &symx); break; } h = (struct ahead *) sym[symy].spec.array.ptr; qybase.l = (int *) ((char *)h + sizeof(struct ahead)); h = (struct ahead *) sym[symx].spec.array.ptr; qxbase.l = (int *) ((char *)h + sizeof(struct ahead)); /*get result sym */ result_sym = array_clone(symy,toptype); h = (struct ahead *) sym[result_sym].spec.array.ptr; qzbase.l = (int *) ((char *)h + sizeof(struct ahead)); /* get the various matrices and vectors needed */ /*unlike pit, the coeffs. are not in a symbol */ if ( (cfbase = (double *) malloc((npow+1) * sizeof(double))) == NULL) return execute_error(23); if ( (fbase = (double *) malloc( nxq * sizeof(double))) == NULL) return execute_error(23); if ( (a = (double *) malloc((npow+1)*(npow+1) * sizeof(double))) == NULL) return execute_error(23); /*loop over the number of individual fits to do */ while (outer--) { /* get the poly fit coefficients, done in cana_pit */ cana_pit( &qxbase, &qybase, nxq, npow, a, cfbase, fbase, toptype, 1); /*cfbase has coeffs., use scratch x to get trend */ cana_poly(cfbase, &qxbase, nxq, npow, &qzbase, toptype); /* check if this is actually a detrend call */ if (detrend_flag) { /* a detrend means subtracting the fit from the data and returning the result rather than the fit */ switch (toptype) { case 3: qybase.f -= nxq; qzbase.f -= nxq; nd = nxq; while (nd--) *qzbase.f++ = *qybase.f++ - *qzbase.f; break; case 4: qybase.d -= nxq; qzbase.d -= nxq; nd = nxq; while (nd--) *qzbase.d++ = *qybase.d++ - *qzbase.d; break; } } } /*end of while (outer--) loop */ free(a); free(fbase); free(cfbase); detrend_flag = 0; return result_sym; } /*------------------------------------------------------------------------- */ int cana_poly( cfbase, qxbase, nxq, npow, qzbase, toptype ) /* internal routine used by ana_trend, and ana_detrend */ union types_ptr *qxbase,*qzbase; double *cfbase; int nxq, npow, toptype; { double sum, *cf, xq; union types_ptr qz,qx; int i, n, k, m; cf = cfbase; qz.l = qzbase->l; qx.l = qxbase->l; switch (toptype) { case 3: n = nxq; switch (npow) { case 0: while (n--) *qz.f++ = *cfbase; break; default: while (n--) { cf = cfbase + npow; /*point to last coeff. */ xq = *qx.f++; sum = (*cf--) * xq; m = npow - 1; while (m-- > 0) sum = ( sum + *cf--) * xq; *qz.f++ = sum + *cf; } } qzbase->f += nxq; break; case 4: n = nxq; switch (npow) { case 0: while (n--) *qz.d++ = *cfbase; break; default: n = nxq; while (n--) { cf = cfbase + npow; /*point to last coeff. */ xq = *qx.d++; sum = (*cf--) * xq; m = npow - 1; while (m-- > 0) sum = ( sum + *cf--) * xq; *qz.d++ = sum + *cf; } } qzbase->d += nxq; break; } return 1; } /*------------------------------------------------------------------------- */ int ana_polate(narg,ps) /*polate function */ /* linear interpolation, CALL IS YF=POLATE(X,Y,XF)*/ int narg,ps[]; { union types_ptr qx,qy,qxf,qyf; double sum, *a, *f, *cf; int npow, symx, symy, symxf, result_sym, nd, nxq, outer, outerx; int iq, jq; int j, i, n, ib, ie, toptype, k, nxf, outerxf; struct ahead *h; /* for 2-D arrays, we do a series of 1-D interpolations */ /* a 2-D y with a 1-D x is permitted, x is used for all the outer y dimensions in this case */ symx = ps[0]; symy = ps[1]; symxf = ps[2]; /*everybody must be an array */ if (sym[symy].class != 4 || sym[symx].class != 4 || sym[symxf].class != 4) return execute_error(66); /* y attributes, y determines outer */ toptype = sym[symy].type < 3 ? 3 : sym[symy].type; h = (struct ahead *) sym[symy].spec.array.ptr; /*array may be upgraded but y dimensions will be same, get now */ nd = h->ndim; nxq = h->dims[0]; outer = 1; if (nd > 1) for(j=1;jdims[j]; /* x attributes, some must match y */ toptype = sym[symx].type > toptype ? sym[symx].type : toptype; h = (struct ahead *) sym[symx].spec.array.ptr; nd = h->ndim; if (h->dims[0] != nxq) return execute_error(77); outerx = 1; if (nd > 1) for(j=1;jdims[j]; if (outerx != 1 && outerx != outer ) return execute_error(77); /* xf attributes */ toptype = sym[symxf].type > toptype ? sym[symxf].type : toptype; h = (struct ahead *) sym[symxf].spec.array.ptr; /*array may be upgraded but y dimensions will be same, get now */ nd = h->ndim; nxf = h->dims[0]; outerxf = 1; if (nd > 1) for(j=1;jdims[j]; if (outerxf != 1 && outerxf != outer ) return execute_error(77); /*now float or double everybody */ switch (toptype) { case 3: symxf = ana_float( 1, &symxf); symy = ana_float( 1, &symy); symx = ana_float( 1, &symx); break; case 4: symxf = ana_double( 1, &symxf); symy = ana_double( 1, &symy); symx = ana_double( 1, &symx); break; } /*get result array */ result_sym = array_clone(symxf, toptype); /*get addresses of everybody */ h = (struct ahead *) sym[symy].spec.array.ptr; qy.l = (int *) ((char *)h + sizeof(struct ahead)); h = (struct ahead *) sym[symxf].spec.array.ptr; qxf.l = (int *) ((char *)h + sizeof(struct ahead)); h = (struct ahead *) sym[symx].spec.array.ptr; qx.l = (int *) ((char *)h + sizeof(struct ahead)); h = (struct ahead *) sym[result_sym].spec.array.ptr; qyf.l = (int *) ((char *)h + sizeof(struct ahead)); /*loop over the number of individual fits to do */ while (outer--) { switch (toptype) { case 3: cana_polate ( nxq, qx.f, qy.f, qxf.f, qyf.f, nxf ); qy.f += nxq; if ( outerx != 1 ) qx.f += outerx; qyf.f += nxq; if ( outerxf != 1 ) qxf.f += outerxf; break; case 4: cana_dpolate ( nxq, qx.d, qy.d, qxf.d, qyf.d, nxf ); qy.d += nxq; if ( outerx != 1 ) qx.d += outerx; qyf.d += nxq; if ( outerxf != 1 ) qxf.d += outerxf; break; } } return result_sym; } /*------------------------------------------------------------------------- */ int cana_polate( nx, qx, qy, qxf, qyf, nxf ) float *qy, *qx, *qxf, *qyf; int nx, nxf; { int n, m, iq, jq, mq, iqold; float dx, grad, xnew; iq = 1; jq = 0; iqold = iq; dx = *(qx + iq) - *qx; if (dx == 0) grad = 0.0; else grad = (*(qy + iq) - *qy) / dx; /*note - this has effect of using first y value if dx = 0 */ /* 4/14/95 changed mq = nx -2 to mg = nx -1; apparent bug */ m = nxf; mq = nx - 1; while (m--) { xnew = *qxf++; /* find the surrounding x's for this value of xf */ /* first find an x > xnew */ /* 12/5/95 yet another bug, the >= below was just > */ while (xnew >= *(qx + iq)) { if (iq >= mq) { iq = mq; break;} iq++; } /* check that we span */ jq = iq -1; while (xnew < *(qx + jq )) { if (jq <= 0 ) { jq = 0; break;} jq--; } iq = jq + 1; if ( iq != iqold ) { dx = *(qx + iq) - *(qx + jq); /* former bug in next line fixed by LS sometime in 9/92 */ if ( dx == 0 ) grad = 0; else grad = (*(qy + iq) - *(qy + jq)) / dx; } *qyf++ = *(qy + jq) + grad * ( xnew - *(qx + jq) ); } return; } /*------------------------------------------------------------------------- */ int cana_dpolate( nx, qx, qy, qxf, qyf, nxf ) double *qy, *qx, *qxf, *qyf; int nx, nxf; { int n, m, iq, jq, mq, iqold; double dx, grad, xnew; iq = 1; jq = 0; iqold = iq; dx = *(qx + iq) - *qx; if (dx == 0) grad = 0.0; else grad = (*(qy + iq) - *qy) / dx; /*note - this has effect of using first y value if dx = 0 */ /* 4/14/95 changed mq = nx -2 to mg = nx -1; apparent bug */ m = nxf; mq = nx - 1; while (m--) { xnew = *qxf++; /* find the surrounding x's for this value of xf */ /* first find an x > xnew */ while (xnew >= *(qx + iq)) { if (iq >= mq) { iq = mq; break;} iq++; } /* check that we span */ jq = iq -1; while (xnew < *(qx + jq )) { if (jq <= 0 ) { jq = 0; break;} jq--; } iq = jq + 1; if ( iq != iqold ) { dx = *(qx + iq) - *(qx + jq); /* former bug in next line fixed by LS sometime in 9/92 */ if ( dx == 0 ) grad = 0; else grad = (*(qy + iq) - *(qy + jq)) / dx; } *qyf++ = *(qy + jq) + grad * ( xnew - *(qx + jq) ); } return; } /*------------------------------------------------------------------------- */ int ana_poly(narg,ps) /*poly function */ /* polynomial generartor, Y = POLY(X,C)*/ int narg,ps[]; { union types_ptr qxbase,qybase,cbase; int npow, symx, symc, nd, nx, outer, outerx, iq, jq, cfsym; int j, i, n, ib, ie, toptype, k, outerc, dim[8]; struct ahead *h; /* if x is 2D or more, the outer dimensions are used as a repeat count (which doesn't matter for this function anyhow); if c is 2D or more, the outers are also assumed to be a repeat count; if both are, then they must match; the result then has the inner dimension of x and the outer of either x or c */ symx = ps[0]; if ( sym[symx].class != 4 ) return execute_error(66); /* x attributes */ toptype = MAX(sym[symx].type, 3); h = (struct ahead *) sym[symx].spec.array.ptr; nd = h->ndim; nx = h->dims[0]; outerx = 1; if (nd > 1) for(j=1;jdims[j]; /* c attributes */ symc = ps[1]; symc = ana_double( 1, &symc ); /*always make double */ outerc = 1; switch (sym[symc].class) { /*scalar c allowed */ case 8: symc = class8_to_1(symc); /*scalar ptr case */ case 1: cbase.l = &sym[symc].spec.scalar.l; npow = 0; break; case 4: h = (struct ahead *) sym[symc].spec.array.ptr; nd = h->ndim; npow = h->dims[0] - 1; if (nd > 1) for(j=1;jdims[j]; cbase.l = (int *) ((char *)h + sizeof(struct ahead)); break; default: return execute_error(83); break; } /* check x and c for consistency and decide on outer for result */ /* x detremines first dimension of result always */ outer = outerx; if (outerc != 1 && outerc != outerx) return execute_error(77); h = (struct ahead *) sym[symx].spec.array.ptr; nd = h->ndim; for(j=0;jdims[j]; if (outerc != 1) { outer = outerc; if (outerx != 1 && outerc != outerx) return execute_error(77); h = (struct ahead *) sym[symc].spec.array.ptr; nd = h->ndim; dim[0] = nx; for(j=1;jdims[j]; } /* upgrade x if necessary */ switch (toptype) { case 3: symx = ana_float( 1, &symx); break; case 4: symx = ana_double( 1, &symx); break; } h = (struct ahead *) sym[symx].spec.array.ptr; qxbase.l = (int *) ((char *)h + sizeof(struct ahead)); /*get result sym */ result_sym = array_scratch(toptype,nd,dim); h = (struct ahead *) sym[result_sym].spec.array.ptr; qybase.l = (int *) ((char *)h + sizeof(struct ahead)); /* do the calculations, looping over outer */ while (outer--) { cana_poly(cbase.d, &qxbase, nx, npow, &qybase, toptype); if (outerc != 1) cbase.d += npow + 1; if (outerx != 1) { switch (toptype) { case 4: qxbase.d += nx; break; case 3: qxbase.f += nx; break; } } } return result_sym; } /*------------------------------------------------------------------------- */ int ana_bilinxy(narg,ps) /*bilinxy function */ /* bi-linear interpolation with a function defined on a rectangular grid a function, result = bilinxy(xim, gx, gy, xf, yf) only floating point supported, no f*8, all I*n types converted to f*4 */ int narg,ps[]; { float *xf, *yf, *out, *xim, *gx, *gy, *pq; float dx, dy, gradx, grady, gxa, gya, xnew, ynew, aq, bq, f00, f01, f10, f11; int iq, jq, kq, iqold, jqold, nx, ny, nq, mq, n, nd, j, result_sym; struct ahead *h; /* xim must be a 2-D array and gx and gy must match the x and y dimensions*/ iq = ps[0]; if ( sym[iq].class != 4 ) return execute_error(66); h = (struct ahead *) sym[iq].spec.array.ptr; if ( h->ndim != 2 ) return execute_error(101); nx = h->dims[0]; ny = h->dims[1]; /* float it, sorry no f*8 support yet */ iq = ana_float(1, &ps[0]); h = (struct ahead *) sym[iq].spec.array.ptr; xim = (float *) ((char *)h + sizeof(struct ahead)); /* gx and gy are the positions for which xim is defined each must be a 1-D array */ iq = ps[1]; if ( sym[iq].class != 4 ) return execute_error(107); h = (struct ahead *) sym[iq].spec.array.ptr; if ( h->ndim != 1 ) return execute_error(107); if (nx != h->dims[0]) return execute_error(77); iq = ana_float(1, &ps[1]); h = (struct ahead *) sym[iq].spec.array.ptr; gx = (float *) ((char *)h + sizeof(struct ahead)); iq = ps[2]; if ( sym[iq].class != 4 ) return execute_error(107); h = (struct ahead *) sym[iq].spec.array.ptr; if ( h->ndim != 1 ) return execute_error(107); if (ny != h->dims[0]) return execute_error(77); iq = ana_float(1, &ps[2]); h = (struct ahead *) sym[iq].spec.array.ptr; gy = (float *) ((char *)h + sizeof(struct ahead)); /* now the xf and yf arrays, check and float */ /* they must be the same length, result has dimensions of xf */ iq = ps[3]; if ( sym[iq].class != 4 ) return execute_error(107); iq = ana_float(1, &ps[3]); h = (struct ahead *) sym[iq].spec.array.ptr; xf = (float *) ((char *)h + sizeof(struct ahead)); nd = h->ndim; n = 1; for (j=0;jdims[j]; /* # of elements for xf */ result_sym = array_clone(iq, sym[iq].type); h = (struct ahead *) sym[result_sym].spec.array.ptr; out = (float *) ((char *)h + sizeof(struct ahead)); iq = ps[4]; if ( sym[iq].class != 4 ) return execute_error(107); iq = ana_float(1, &ps[4]); h = (struct ahead *) sym[iq].spec.array.ptr; yf = (float *) ((char *)h + sizeof(struct ahead)); nd = h->ndim; nq = 1; for (j=0;jdims[j]; /* # of elements for yf */ if ( nq != n) return execute_error(77); /* initial dx and dy and iq and jq */ iq = iqold =1; jq = jqold = 1; dx = *(gx + iq) - *gx; dy = *(gy + jq) - *gy; if (dx == 0) gradx = 0.0; else gradx = 1.0 / dx; if (dy == 0) grady = 0.0; else grady = 1.0 / dy; nq = nx - 1; mq = ny - 1; gxa = *gx; gya = *gy; while (n--) { /* step through (xf, yf) and compute out from xim, find bounds on x */ xnew = *xf++; /* find the surrounding x's for this value of xf */ /* first find an x > xnew */ while (xnew >= *(gx + iq)) { if (iq >= nq) { iq = nq; break;} iq++; } /* but may have to go backwards */ kq = iq - 1; while (xnew < *(gx + kq )) { if (kq <= 0 ) { kq = 0; break;} kq--; } iq = kq + 1; if ( iq != iqold ) { gxa = *(gx + kq); dx = *(gx + iq) - gxa; iqold = iq; if ( dx == 0 ) gradx = 0; else gradx = 1.0 / dx; } aq = (xnew - gxa) * gradx; /* now find the surrounds for y */ ynew = *yf++; while (ynew >= *(gy + jq)) { if (jq >= mq) { jq = mq; break;} jq++; } /* but may have to go backwards */ kq = jq - 1; while (ynew < *(gy + kq )) { if (kq <= 0 ) { kq = 0; break;} kq--; } jq = kq + 1; if ( jq != jqold ) { gya = *(gy + kq); dy = *(gy + jq) - gya; jqold = jq; if ( dy == 0 ) grady = 0; else grady = 1.0 / dy; } bq = (ynew - gya) * grady; /* get upper corner index */ /* printf("xnew, ynew = %f %f iq, jq = %d %d\n", xnew, ynew, iq,jq); printf("aq, bq = %f %f\n", aq, bq); */ pq = xim + iq + nx * jq; f11 = *pq; pq--; f10 = *pq; pq = pq -nx; f00 = *pq; pq++; f01 = *pq; /* printf("f00, f01, f10, f11 = %f %f %f %f\n", f00, f01, f10, f11); */ *out++ = aq*bq*(f00+f11-f01-f10) + aq*(f01-f00) + bq*(f10-f00) + f00; } return result_sym; } /*end of ana_bilinxy */ /*------------------------------------------------------------------------- */ int ana_bilinxy_dp(narg,ps) /* bilinxy function, double precision version */ /* bi-linear interpolation with a function defined on a rectangular grid a function, result = bilinxy(xim, gx, gy, xf, yf) all I*n types converted to f*8 */ int narg,ps[]; { double *xf, *yf, *out, *xim, *gx, *gy, *pq; double dx, dy, gradx, grady, gxa, gya, xnew, ynew, aq, bq, f00, f01, f10, f11; int iq, jq, kq, iqold, jqold, nx, ny, nq, mq, n, nd, j, result_sym; struct ahead *h; /* xim must be a 2-D array and gx and gy must match the x and y dimensions*/ iq = ps[0]; if ( sym[iq].class != 4 ) return execute_error(66); h = (struct ahead *) sym[iq].spec.array.ptr; if ( h->ndim != 2 ) return execute_error(101); nx = h->dims[0]; ny = h->dims[1]; /* double it */ iq = ana_double(1, &ps[0]); h = (struct ahead *) sym[iq].spec.array.ptr; xim = (double *) ((char *)h + sizeof(struct ahead)); /* gx and gy are the positions for which xim is defined each must be a 1-D array */ iq = ps[1]; if ( sym[iq].class != 4 ) return execute_error(107); h = (struct ahead *) sym[iq].spec.array.ptr; if ( h->ndim != 1 ) return execute_error(107); if (nx != h->dims[0]) return execute_error(77); iq = ana_double(1, &ps[1]); h = (struct ahead *) sym[iq].spec.array.ptr; gx = (double *) ((char *)h + sizeof(struct ahead)); iq = ps[2]; if ( sym[iq].class != 4 ) return execute_error(107); h = (struct ahead *) sym[iq].spec.array.ptr; if ( h->ndim != 1 ) return execute_error(107); if (ny != h->dims[0]) return execute_error(77); iq = ana_double(1, &ps[2]); h = (struct ahead *) sym[iq].spec.array.ptr; gy = (double *) ((char *)h + sizeof(struct ahead)); /* now the xf and yf arrays, check and float */ /* they must be the same length, result has dimensions of xf */ iq = ps[3]; if ( sym[iq].class != 4 ) return execute_error(107); iq = ana_double(1, &ps[3]); h = (struct ahead *) sym[iq].spec.array.ptr; xf = (double *) ((char *)h + sizeof(struct ahead)); nd = h->ndim; n = 1; for (j=0;jdims[j]; /* # of elements for xf */ result_sym = array_clone(iq, sym[iq].type); h = (struct ahead *) sym[result_sym].spec.array.ptr; out = (double *) ((char *)h + sizeof(struct ahead)); iq = ps[4]; if ( sym[iq].class != 4 ) return execute_error(107); iq = ana_double(1, &ps[4]); h = (struct ahead *) sym[iq].spec.array.ptr; yf = (double *) ((char *)h + sizeof(struct ahead)); nd = h->ndim; nq = 1; for (j=0;jdims[j]; /* # of elements for yf */ if ( nq != n) return execute_error(77); /* initial dx and dy and iq and jq */ iq = iqold =1; jq = jqold = 1; dx = *(gx + iq) - *gx; dy = *(gy + jq) - *gy; if (dx == 0) gradx = 0.0; else gradx = 1.0 / dx; if (dy == 0) grady = 0.0; else grady = 1.0 / dy; nq = nx - 1; mq = ny - 1; gxa = *gx; gya = *gy; while (n--) { /* step through (xf, yf) and compute out from xim, find bounds on x */ xnew = *xf++; /* find the surrounding x's for this value of xf */ /* first find an x > xnew */ while (xnew >= *(gx + iq)) { if (iq >= nq) { iq = nq; break;} iq++; } /* but may have to go backwards */ kq = iq - 1; while (xnew < *(gx + kq )) { if (kq <= 0 ) { kq = 0; break;} kq--; } iq = kq + 1; if ( iq != iqold ) { gxa = *(gx + kq); dx = *(gx + iq) - gxa; iqold = iq; if ( dx == 0 ) gradx = 0; else gradx = 1.0 / dx; } aq = (xnew - gxa) * gradx; /* now find the surrounds for y */ ynew = *yf++; while (ynew >= *(gy + jq)) { if (jq >= mq) { jq = mq; break;} jq++; } /* but may have to go backwards */ kq = jq - 1; while (ynew < *(gy + kq )) { if (kq <= 0 ) { kq = 0; break;} kq--; } jq = kq + 1; if ( jq != jqold ) { gya = *(gy + kq); dy = *(gy + jq) - gya; jqold = jq; if ( dy == 0 ) grady = 0; else grady = 1.0 / dy; } bq = (ynew - gya) * grady; /* get upper corner index */ /* printf("xnew, ynew = %f %f iq, jq = %d %d\n", xnew, ynew, iq,jq); printf("aq, bq = %f %f\n", aq, bq); */ pq = xim + iq + nx * jq; f11 = *pq; pq--; f10 = *pq; pq = pq -nx; f00 = *pq; pq++; f01 = *pq; /* printf("f00, f01, f10, f11 = %f %f %f %f\n", f00, f01, f10, f11); */ *out++ = aq*bq*(f00+f11-f01-f10) + aq*(f01-f00) + bq*(f10-f00) + f00; } return result_sym; } /*end of ana_bilinxy */ /*----------------------------------------------------------------------- */