From fa080de7afc95aa1c19a6e6fc0e0708ced2eadc4 Mon Sep 17 00:00:00 2001 From: Joseph Hunkeler Date: Wed, 8 Jul 2015 20:46:52 -0400 Subject: Initial commit --- noao/obsutil/src/sptime/tabinterp.x | 698 ++++++++++++++++++++++++++++++++++++ 1 file changed, 698 insertions(+) create mode 100644 noao/obsutil/src/sptime/tabinterp.x (limited to 'noao/obsutil/src/sptime/tabinterp.x') diff --git a/noao/obsutil/src/sptime/tabinterp.x b/noao/obsutil/src/sptime/tabinterp.x new file mode 100644 index 00000000..518f7b20 --- /dev/null +++ b/noao/obsutil/src/sptime/tabinterp.x @@ -0,0 +1,698 @@ +include +include +include + +# Table structure. +define TAB_DIM 2 # Maximum dimension of table +define TAB_SZFNAME 99 # Size of file name +define TAB_SZPARAMS (10*SZ_LINE) # Size of parameter string +define TAB_SIZE (55+2*TAB_DIM) # Size of table structure + +define TAB_FNAME Memc[P2C($1)] # File name +define TAB_VALUE Memr[P2R($1+50)] # Constant table value +define TAB_PARAMS Memi[$1+51] # Pointer to parameters +define TAB_NDIM Memi[$1+52] # Dimension of table +define TAB_INTERP Memi[$1+53] # Interpolation pointer +define TAB_NEXTRAP Memi[$1+54] # Number of extrapolations +define TAB_LEN Memi[$1+54+$2] # Length of axes in table +define TAB_COORD Memi[$1+54+$2+TAB_DIM] # Pointer to axes coordinates + + +procedure tabtest () + +char table[SZ_FNAME], param[SZ_FNAME], str[SZ_FNAME] +int clglpr() +real x, y, z, tabinterp1(), tabinterp2(), tabgetr() +pointer tab, tabopen() +errchk tabopen, tabinterp1, tabinterp2, tabgetr, tabgstr + +begin + tab = tabopen () + + call clgstr ("table", table, SZ_FNAME) + + call clgstr ("param", param, SZ_FNAME) + z = tabgetr (tab, table, "", param) + call tabgstr (tab, table, "", param, str, SZ_FNAME) + call printf ("%s = %g = %s\n") + call pargstr (param) + call pargr (z) + call pargstr (str) + + while (clglpr ("x", x) != EOF) { + iferr { + z = tabinterp1 (tab, table, x) + call printf ("%g %g\n") + call pargr (x) + call pargr (z) + } then { + if (clglpr ("y", y) == EOF) + break + z = tabinterp2 (tab, table, x, y) + call printf ("%g %g %g\n") + call pargr (x) + call pargr (y) + call pargr (z) + } + } + + call tabclose (tab) +end + + +# TABOPEN -- Open table interpolation package. + +pointer procedure tabopen () + +pointer stp, stopen() + +begin + stp = stopen ("tables", 10, 1000, 1000) + return (stp) +end + + +# TABCLOSE -- Close table interpolation package. + +procedure tabclose (stp) + +pointer stp #I Symbol table pointer + +pointer sym, sthead(), stnext() + +begin + for (sym = sthead(stp); sym != NULL; sym = stnext(stp, sym)) + call tabfree (sym) + call stclose (stp) +end + + +procedure tabfree (sym) + +pointer sym #I Table structure + +int i + +begin + call mfree (TAB_PARAMS(sym), TY_CHAR) + if (TAB_INTERP(sym) != NULL) { + if (TAB_LEN(sym,1) > 1 && TAB_LEN(sym,2) > 1) + call msifree (TAB_INTERP(sym)) + else + call asifree (TAB_INTERP(sym)) + do i = 1, TAB_NDIM(sym) + call mfree (TAB_COORD(sym,i), TY_REAL) + } +end + + + +# TABLOAD -- Load a table in the symbol table. + +procedure tabload (stp, name, fname) + +pointer stp # Symbol table pointer +char name[ARB] # Name of table +char fname[ARB] # File name of table + +int i, fd, ndim, npts, len[2], open(), errget(), nowhite(), ctor(), strlen() +real value +pointer sp, str, sym, params, data, coords[2], stfind(), stenter() +bool streq() +errchk open, tabread + +begin + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + + # If no value is specified then don't enter into symbol table. + if (nowhite (fname, Memc[str], SZ_LINE) == 0) { + call sfree (sp) + return + } + + # The special string "none" is equivalent to no table. + if (streq (Memc[str], "none")) { + call sfree (sp) + return + } + + # Check if table has already been loaded. + sym = stfind (stp, name) + if (sym != NULL) { + if (streq (fname, TAB_FNAME(sym))) + return + } + + # Check if constant is specified. + i = 1 + if (ctor (Memc[str], i, value) == strlen (Memc[str])) { + sym = stenter (stp, name, TAB_SIZE) + call strcpy (fname, TAB_FNAME(sym), TAB_SZFNAME) + TAB_VALUE(sym) = value + TAB_NDIM(sym) = -1 + TAB_PARAMS(sym) = NULL + TAB_INTERP(sym) = NULL + + call sfree (sp) + return + } + + # Read the table. + fd = open (Memc[str], READ_ONLY, TEXT_FILE) + iferr (call tabread (fd, params, data, npts, len, coords, ndim)) { + ndim = errget (Memc[str], SZ_LINE) + call strcat (" (", Memc[str], SZ_LINE) + call strcat (name, Memc[str], SZ_LINE) + call strcat (")", Memc[str], SZ_LINE) + call error (1, Memc[str]) + } + call close (fd) + + if (sym == NULL) + sym = stenter (stp, name, TAB_SIZE) + else + call tabfree (sym) + + if (data != NULL) { + if (len[1] > 1 && len[2] > 1) { + call msiinit (TAB_INTERP(sym), II_BILINEAR) + call msifit (TAB_INTERP(sym), Memr[data], len[1], len[2], + len[1]) + } else if (len[2] == 1) { + if (len[1] == 1) + call asiinit (TAB_INTERP(sym), II_NEAREST) + else + call asiinit (TAB_INTERP(sym), II_LINEAR) + call asifit (TAB_INTERP(sym), Memr[data], len[1]) + } else { + call asiinit (TAB_INTERP(sym), II_LINEAR) + call asifit (TAB_INTERP(sym), Memr[data], len[2]) + } + } else + TAB_INTERP(sym) = NULL + + call strcpy (fname, TAB_FNAME(sym), TAB_SZFNAME) + TAB_PARAMS(sym) = params + TAB_NDIM(sym) = ndim + TAB_NEXTRAP(sym) = 0 + call amovi (len, TAB_LEN(sym,1), 2) + call amovi (coords, TAB_COORD(sym,1), 2) + call mfree (data, TY_REAL) + call sfree (sp) +end + + +# TABREAD -- Read data from a table. + +procedure tabread (fd, params, data, npts, len, coords, ndim) + +int fd #I File descriptor +pointer params #O Pointer to parameters +pointer data #O Pointer to data +int npts #O Number of data points +int len[ARB] #O Number of points along each dimension +pointer coords[ARB] #O Pointers to coordinates along each dimension +int ndim #O Dimension of table + +int i, j, nread, fdparams, fscan(), nscan(), stropen(), ctor() +real coord[4], scale +pointer sp, cmt, key, eq, val +bool streq() +errchk stropen() + +begin + iferr { + call smark (sp) + call salloc (cmt, 2, TY_CHAR) + call salloc (key, SZ_FNAME, TY_CHAR) + call salloc (eq, 1, TY_CHAR) + call salloc (val, SZ_LINE, TY_CHAR) + + npts = 0 + ndim = 0 + params = NULL + data = NULL + scale = INDEF + call aclri (len, 2) + call aclri (coords, 2) + while (fscan (fd) != EOF) { + call gargr (coord[1]) + call gargr (coord[2]) + call gargr (coord[3]) + call gargr (coord[4]) + nread = nscan() + + if (nread == 0) { # Check for parameters. + call reset_scan() + call gargwrd (Memc[cmt], 2) + call gargwrd (Memc[key], SZ_LINE) + call gargwrd (Memc[eq], 1) + call gargwrd (Memc[val], SZ_LINE) + if (nscan() != 4 || Memc[cmt] != '#' || Memc[eq] != '=' || + Memc[cmt+1] != EOS) + next + if (streq (Memc[key], "tablescale")) { + i = 1 + if (ctor (Memc[val], i, scale) == 0) + call error (1, "Syntax error in table scale factor") + } + if (params == NULL) { + call malloc (params, TAB_SZPARAMS, TY_CHAR) + fdparams = stropen (Memc[params], TAB_SZPARAMS, + WRITE_ONLY) + } + call fprintf (fdparams, "%s %s\n") + call pargstr (Memc[key]) + call pargstr (Memc[val]) + next + } else if (nread == 1) + next + + if (ndim == 0) { + ndim = nread - 1 + if (ndim > 2) + call error (1, "Table dimension is too high") + do i = ndim+1, 2 + len[i] = 1 + } + if (nread-1 != ndim) + call error (2, "Table has variable number of columns") + + do i = 1, ndim { + if (len[i] == 0) { + call malloc (coords[i], 100, TY_REAL) + Memr[coords[i]] = coord[i] + len[i] = len[i] + 1 + } else { + do j = 0, len[i] - 1 + if (coord[i] == Memr[coords[i]+j]) + break + if (j >= len[i]) { + if (mod (len[i], 100) == 0) + call realloc (coords[i], len[i]+100, TY_REAL) + Memr[coords[i]+len[i]] = coord[i] + len[i] = len[i] + 1 + } + } + } + + if (npts == 0) + call malloc (data, 100, TY_REAL) + else if (mod (npts, 100) == 0) + call realloc (data, npts+100, TY_REAL) + + Memr[data+npts] = coord[nread] + npts = npts + 1 + } + + if (npts > 0) { + j = 1 + do i = 1, ndim + j = j * len[i] + if (j != npts) + call error (4, "Table is not regular") + } + + if (!IS_INDEF(scale)) + call amulkr (Memr[data], scale, Memr[data], npts) + + call close (fdparams) + call sfree (sp) + } then { + if (params != NULL) { + call close (fdparams) + call mfree (params, TY_CHAR) + } + do i = 1, 2 + call mfree (coords[i], TY_REAL) + call mfree (data, TY_REAL) + call sfree (sp) + call erract (EA_ERROR) + } +end + + +# TABEXISTS -- Determine if table exists. + +bool procedure tabexists (stp, name) + +pointer stp #I Symbol table pointer +char name[ARB] #I Name of table + +pointer stfind() + +begin + return (stfind (stp, name) != NULL) +end + + +# TABGETR -- Get real parameter from table. + +real procedure tabgetr (stp, name, alt, param) + +pointer stp #I Symbol table pointer +char name[ARB] #I Name of table +char alt[ARB] #I Name of alternate table +char param[ARB] #I Parameter name +real val #O Returned value + +real rval +int i, fd, strncmp(), strdic(), stropen(), fscan(), nscan() +bool streq() +pointer sp, key, sym[2], stfind() +errchk stfind, stropen + +begin + # Return if no table. + if (name[1] == EOS && alt[1] == EOS) + return (INDEFR) + + call smark (sp) + call salloc (key, SZ_FNAME, TY_CHAR) + + # Find tables. + sym[1] = stfind (stp, name) + if (alt[1] != EOS) + sym[2] = stfind (stp, alt) + else + sym[2] = NULL + + # Return an error if there is no table. + if (sym[1] == NULL && sym[2] == NULL) { + call sprintf (Memc[key], SZ_FNAME, "Table `%s' not found") + call pargstr (name) + call error (1, Memc[key]) + } + + # Get the parameter value. + val = INDEFR + do i = 1, 2 { + if (sym[i] == NULL) + next + if (strncmp (param, "table.", 6) == 0) { + switch (strdic (param[7], Memc[key], SZ_FNAME, + "|ndim|xmin|xmax|ymin|ymax|nextrap|")) { + case 1: + val = TAB_NDIM(sym[i]) + case 2: + if (TAB_NDIM(sym[i]) >= 1) + val = Memr[TAB_COORD(sym[i],1)] + else if (TAB_NDIM(sym[i]) == -1) + val = -MAX_REAL + case 3: + if (TAB_NDIM(sym[i]) >= 1) + val = Memr[TAB_COORD(sym[i],1)+TAB_LEN(sym[i],1)-1] + else if (TAB_NDIM(sym[i]) == -1) + val = MAX_REAL + case 4: + if (TAB_NDIM(sym[i]) >= 2) + val = Memr[TAB_COORD(sym[i],2)] + else if (TAB_NDIM(sym[i]) == -1) + val = -MAX_REAL + case 5: + if (TAB_NDIM(sym[i]) >= 2) + val = Memr[TAB_COORD(sym[i],2)+TAB_LEN(sym[i],1)-1] + else if (TAB_NDIM(sym[i]) == -1) + val = MAX_REAL + case 6: + val = TAB_NEXTRAP(sym[i]) + } + break + + } else if (TAB_PARAMS(sym[i]) != NULL) { + fd = stropen (Memc[TAB_PARAMS(sym[i])], TAB_SZPARAMS, READ_ONLY) + while (fscan (fd) != EOF) { + call gargwrd (Memc[key], SZ_FNAME) + call gargr (rval) + if (nscan() != 2) + next + if (streq (Memc[key], param)) { + val = rval + break + } + } + call close (fd) + } + + if (!IS_INDEF(val)) + break + } + + # Return error if no value was found. + if (IS_INDEF(val)) { + call sprintf (Memc[key], SZ_FNAME, + "Table parameter `%s' not found (%s)") + call pargstr (param) + call pargstr (name) + call error (1, Memc[key]) + } + + call sfree (sp) + return (val) + +end + + +# TABGETI -- Get integer paraemter from table. + +int procedure tabgeti (stp, name, alt, param) + +pointer stp #I Symbol table pointer +char name[ARB] #I Name of table +char alt[ARB] #I Name of alternate table +char param[ARB] #I Parameter name + +real val, tabgetr() +errchk tabgetr + +begin + val = tabgetr (stp, name, alt, param) + if (IS_INDEFR(val)) + return (INDEFI) + + return (nint(val)) +end + + +# TABGSTR -- Get string parameter from table. + +procedure tabgstr (stp, name, alt, param, val, maxchar) + +pointer stp #I Symbol table pointer +char name[ARB] #I Name of table +char alt[ARB] #I Name of alternate table +char param[ARB] #I Parameter name +char val[ARB] #O Returned value +int maxchar #I Maximum string length + +int i, fd, strncmp(), strdic(), stropen(), fscan(), nscan() +bool streq() +pointer sp, key, str, sym[2], stfind() +errchk stfind, stropen() + +begin + # Return if no table. + if (name[1] == EOS && alt[1] == EOS) { + val[1] = EOS + return + } + + call smark (sp) + call salloc (key, SZ_FNAME, TY_CHAR) + call salloc (str, maxchar, TY_CHAR) + + # Find tables. + sym[1] = stfind (stp, name) + if (alt[1] != EOS) + sym[2] = stfind (stp, alt) + else + sym[2] = NULL + + # Return an error if there is no table. + if (sym[1] == NULL && sym[2] == NULL) { + call sprintf (Memc[key], SZ_FNAME, "Table `%s' not found") + call pargstr (name) + call error (1, Memc[key]) + } + + # Get the parameter value. + val[1] = EOS + do i = 1, 2 { + if (sym[i] == NULL) + next + if (strncmp (param, "table.", 6) == 0) { + switch (strdic (param[7], Memc[key], SZ_FNAME, "|filename|")) { + case 1: + call strcpy (TAB_FNAME(sym[i]), val, maxchar) + } + break + + } else if (TAB_PARAMS(sym[i]) != NULL) { + fd = stropen (Memc[TAB_PARAMS(sym[i])], TAB_SZPARAMS, READ_ONLY) + while (fscan (fd) != EOF) { + call gargwrd (Memc[key], SZ_FNAME) + call gargstr (Memc[str], SZ_LINE) + if (nscan() != 2) + next + if (streq (Memc[key], param)) { + if (val[1] == EOS) + call strcpy (Memc[str+1], val, maxchar) + else { + call strcat ("\n", val, maxchar) + call strcat (Memc[str+1], val, maxchar) + } + } + } + call close (fd) + + if (val[1] != EOS) + break + } + } + + # Return error if no value was found. + if (val[1] == EOS) { + call sprintf (Memc[key], SZ_FNAME, + "Table parameter `%s' not found (%s)") + call pargstr (param) + call pargstr (name) + call error (1, Memc[key]) + } + + call sfree (sp) +end + + +# TABINTERP1 -- Interpolate a named 1D table. + +real procedure tabinterp1 (stp, name, x) + +pointer stp # Symbol table pointer +char name[ARB] # Name of table +real x # Interpolation coordinate + +char err[SZ_FNAME] +int i, nx +real xi, y, asieval() +pointer sym, xp, stfind() +errchk stfind + +begin + # Find the table. + sym = stfind (stp, name) + if (sym == NULL) { + call sprintf (err, SZ_FNAME, "Table `%s' not found") + call pargstr (name) + call error (1, err) + } + + # If a constant table return the value. + if (TAB_NDIM(sym) == -1) + return (TAB_VALUE(sym)) + + # Check if the table is of the proper dimensionality. + if (TAB_NDIM(sym) != 1) { + call sprintf (err, SZ_FNAME, "Table is not one dimensional (%s)") + call pargstr (name) + call error (1, err) + } + + nx = TAB_LEN(sym,1) + xp = TAB_COORD(sym,1) + if (x < Memr[xp] || x > Memr[xp+nx-1]) + TAB_NEXTRAP(sym) = TAB_NEXTRAP(sym) + 1 + + if (nx == 1) + xi = 1 + else { + do i = 1, nx-2 + if (x < Memr[xp+i]) + break + + xi = i + (x - Memr[xp+i-1]) / (Memr[xp+i] - Memr[xp+i-1]) + } + xi = max (1., min (real(nx), xi)) + y = asieval (TAB_INTERP(sym), xi) + + + return (y) +end + + +# TABINTERP2 -- Interpolate a named 2D table. + +real procedure tabinterp2 (stp, name, x, y) + +pointer stp # Symbol table pointer +char name[ARB] # Name of table +real x, y # Interpolation coordinate + +char err[SZ_FNAME] +int i, nx, ny +real xi, yi, z, asieval(), msieval() +pointer sym, xp, yp, stfind() +errchk stfind + +begin + # Find the table. + sym = stfind (stp, name) + if (sym == NULL) { + call sprintf (err, SZ_FNAME, "Table `%s' not found") + call pargstr (name) + call error (1, err) + } + + # If a constant table return the value. + if (TAB_NDIM(sym) == -1) + return (TAB_VALUE(sym)) + + # Check if the table is of the proper dimensionality. + if (TAB_NDIM(sym) != 2) { + call sprintf (err, SZ_FNAME, "Table is not two dimensional (%s)") + call pargstr (name) + call error (1, err) + } + + nx = TAB_LEN(sym,1) + ny = TAB_LEN(sym,2) + if (nx > 1 && ny > 1) { # 2D interpolation + xp = TAB_COORD(sym,1) + do i = 1, nx-2 + if (x < Memr[xp+i]) + break + xi = i + (x - Memr[xp+i-1]) / (Memr[xp+i] - Memr[xp+i-1]) + xi = max (1., min (real(nx), xi)) + yp = TAB_COORD(sym,2) + do i = 1, ny-2 + if (y < Memr[yp+i]) + break + yi = i + (y - Memr[yp+i-1]) / (Memr[yp+i] - Memr[yp+i-1]) + yi = max (1., min (real(nx), yi)) + z = msieval (TAB_INTERP(sym), xi, yi) + } else if (ny == 1) { # 1D interpolation in x + if (nx == 1) + xi = 1 + else { + xp = TAB_COORD(sym,1) + do i = 1, nx-2 + if (x < Memr[xp+i]) + break + + xi = i + (x - Memr[xp+i-1]) / (Memr[xp+i] - Memr[xp+i-1]) + } + xi = max (1., min (real(nx), xi)) + z = asieval (TAB_INTERP(sym), xi) + } else { # 1D interpolation in y + yp = TAB_COORD(sym,2) + do i = 1, ny-2 + if (y < Memr[yp+i]) + break + + yi = i + (y - Memr[yp+i-1]) / (Memr[yp+i] - Memr[yp+i-1]) + yi = max (1., min (real(ny), yi)) + z = asieval (TAB_INTERP(sym), yi) + } + + return (z) +end -- cgit