local yodel; /* DOCUMENT yodel.i --- IDL-compatible functions List of routines * strnum: convert number to a "nice" string * strrepl: replace character occurrences within string * strpad: pad string * strpos: locate substring within a string * strmid: extract substring * strsplit: split string using delimiter * strcombine: combine strings using delimiter * strloc: locate given string in a list of strings * rangeloc: locate numeric/string range in within a list of values * dim_conform: check if two dimension lists conform * dim_reshape: reshape dimension list * array_eq: check if two arrays are equal * reshape_array: reshape array dimensions * alleq: check if all elements of an array are equal to each other * monotonic: check if array elements are strictly monotonic * broadcast: broadcast array dimensions * rangeop: carry out operations on selected dimension * arrayop: efficient version of RANGEOP for 5D arrays (for certain ops) * xrotate: periodic coordinate rotation operator * lonlabel: longitude label generator * latlabel: latitude label generator * datelabel: date label generator * param_set: check if parameter is set to a non-null/non-zero value * imaginary: return imaginary part of complex number * is_null: check if object is null * is_integer: check if object is of integer type * is_number: check if object is of a numeric type * is_where: check if "where" output is non-null * prod: return product of all elements of an array * outer_prod: compute outer product of two vectors * low_pass_filt: compute weights for a nonrecursive low-pass filter * randomu1: return uniformly-distributed random number ove unit interval * ref: return reference to ("address of") non-null object * deref: dereference non-null pointer * typeconv: convert object to specified type * raw_dump: read/write raw binary data * oscommand: execute operating system command * timer_init: execution time counter (initialization) * timer_call: execution time counter (function call) * timer_return: execution time counter (function return) * timer_show: execution time counter (display) * * Version: 1.3 * Modified: 2 Mar 1999, R. Saravanan * * SEE ALSO: flexplot */ require, "string.i" func strnum(number,format=) /* DOCUMENT strnum(number,format=) Converts number to "nice" looking string, with optional formatting. SEE ALSO: strrepl, strcombine, strpos, strmid */ { func_name= "strnum"; timer_call,func_name; if (dimsof(number)(1) > 0) { // Array of numbers; handle recursively str= array(string, dimsof(number)) for (i=1; i <= numberof(number); i++) { str(i)= strnum(number(i),format=format); } return timer_return(func_name,str); } // Convert number to string if (is_void(format)) str= swrite(number); else str= swrite(number,format=format); // Remove all spaces, and convert to lower case str= strtolower( strrepl(str,all=1) ); // Extract exponent (if any) expstr= ""; if ( (iexp = strpos(strtolower(str),"e")) != -1) { expstr= strpart(str, iexp+1:); str= strpart(str, :iexp); } if (strpos(str,".") != -1) { // Number contains decimal point; removing trailing zeros str= strrepl(str,'0'); // Strip any trailing decimal point str= strrepl(str,'.'); } // Return (with exponent, if any) return timer_return(func_name,str+expstr); } func strrepl(str,ch,repch,all=) /* DOCUMENT strrepl(str,ch,repch) replace all trailing occurrences of character CH (which defaults to space) in string STR with character REPCH (which defaults to the null character). If ALL==1, then all occurrences of CH are replaced. SEE ALSO: strpad, strnum, strpos, strmid, strsplit */ { if (is_void(ch)) ch= ' '; if (is_void(repch)) repch= '\0'; if (is_void(all)) all= 0; carr= *(pointer(str)); if (all) { if (repch == '\0') carr= carr(where(carr != ch)); else carr(where(carr == ch))= repch; return string(&carr); } else { // Locate trailing occurrence of character len= strlen(str); for (i=len; (i>0) && (carr(i) == ch) ; i--); if (repch != '\0') { if (i < len) carr(i+1:)= repch; return string(&carr); } else { if (i > 0) return strpart(str,1:i); else return ""; } } } func strpad( str, nch, padch, truncate=) /* DOCUMENT strpad(str,nch,padch,truncate=0/1) * Pads string STR with PADCH (which defaults to space) until it is * NCH characters long. * If TRUNCATE==1, if STR is longer than NCH, it is truncated. * SEE ALSO: strrepl, strcombine */ { // Compatibility parameters (including error handling) // compatible.pro: // Parameters for IDL/Yorick compatibility (including error handling) // Array starting index offset, and final index negative offset // (0, 1 in IDL; 1, 0 in Yorick) I0=1 ; I1=0; // Null value ("" in IDL; [] in Yorick) NULL= []; if (!is_scalar(str)) { // Array of strings; handle recursively outstr= str; nstr= numberof(outstr); for (j=I0; j <= nstr-I1; j++) { outstr(j)= strpad(str(j), nch, padch, truncate=truncate); } return outstr; } nlen= strlen(str); if (nlen == nch) return str; if (nlen > nch) { if (param_set(truncate)) return strmid(str,0,nch) ; else return str; } padch1= " "; if (!is_void(padch)) padch1= padch; return str + strcombine(array(padch1,nch-nlen),""); } func strpos(str,pat) /* DOCUMENT strpos(str,pat) locates pattern PAT in string STR, returning -1 if not found, or the starting offset in the string (>=0) where the pattern is found. (A null pattern always returns -1.) SEE ALSO: strnum, strrepl */ { slen= strlen(str); plen= strlen(pat); if ((plen == 0) || (plen > slen)) return -1; if (!strmatch(str,pat)) return -1; for (i= 1; i<=slen-plen+1; i++) if (strpart( str, i:i+plen-1 ) == pat) return i-1; } func strmid(str,offset,length) /* DOCUMENT strmid(str,offset,length) returns the substring of string STR, starting at OFFSET (>=0) with length LENGTH. If there are not enough characters in STR, then the substring is truncated. SEE ALSO: strnum, strrepl */ { if (dimsof(str)(1) > 0) { // Array of strings; handle recursively outstr= array(string, dimsof(str)) for (i=1; i <= numberof(str); i++) { outstr(i)= strmid(str(i),offset,length); } return outstr; } minlen= min([strlen(str)-offset, length]); if (minlen <= 0) return ""; return strpart(str, offset+1:offset+length); } func strsplit( str, delim) /* DOCUMENT strsplit(str, delim) * Decompose string STR into an array of component strings separated by * character DELIM, and return an array of the component strings. * If DELIM is the null string, split string into its constituent one-character * strings. * (Reverses the action of STRCOMBINE) * SEE ALSO: strrepl, strcombine */ { func_name= "strsplit"; timer_call,func_name; // Compatibility parameters (including error handling) // compatible.pro: // Parameters for IDL/Yorick compatibility (including error handling) // Array starting index offset, and final index negative offset // (0, 1 in IDL; 1, 0 in Yorick) I0=1 ; I1=0; // Null value ("" in IDL; [] in Yorick) NULL= []; if (delim == "") { // Split string into individual character strings nlen= strlen(str); if (nlen == 0) return timer_return(func_name, str); str2= array("",nlen); for (j=I0; j <= nlen-I1; j++) str2(j)= strmid( str, j-I0, 1); return timer_return(func_name, str2); } str2= str; str_array= NULL; while ( ((i = strpos(str2,delim))) >= 0 ) { if (i == 0) { grow, str_array, ""; } else { grow, str_array, strmid(str2,0,i); } str2= strmid(str2,i+1,strlen(str2)-i-1); } grow, str_array, str2; return timer_return(func_name, str_array); } func strcombine( str_array, delim) /* DOCUMENT strcombine(str_array,delim) * Concatenate an array of strings STR_ARRAY into a single string, * with each substring separated by character DELIM (which should not occur in * any of the strings. * If DELIM is the null string, simply all concatenate all the strings in * STR_ARRAY. * (Reverses the action of STRSPLIT) * SEE ALSO: strrepl, strsplit */ { func_name= "strcombine"; timer_call,func_name; // Compatibility parameters (including error handling) // compatible.pro: // Parameters for IDL/Yorick compatibility (including error handling) // Array starting index offset, and final index negative offset // (0, 1 in IDL; 1, 0 in Yorick) I0=1 ; I1=0; // Null value ("" in IDL; [] in Yorick) NULL= []; catstr= ""; nstr= numberof(str_array); for (j=I0; j <= nstr-I1; j++) { if (delim != "") { if (strpos(str_array(j), delim) >= 0) error, "Error - delimiter occurs in string"; } if (j == I0) { catstr= str_array(j); } else { catstr= catstr + delim + str_array(j); } } return timer_return(func_name, catstr); } func strloc(name_list,label,case_fold=, abbrev=,comment=) /* DOCUMENT strloc,name_list,label,case_fold=case_fold, * abbrev=abbrev,comment=comment * Locate string LABEL in a list of strings NAME_LIST * return label_index (>=1), if found; else return 0. * If NAME_LIST is null, always return 0. * If LABEL is a null string, always return 0. * If LABEL is an integer index, simply return it, if it is within range. * If CASE_FOLD is true, the matching is case-insensitive * If ABBREV is true, LABEL may be an unambiguous abbreviation * of one of the strings in NAME_LIST. * If COMMENT is non-null, display an error message. * SEE ALSO: strpos, rangeloc */ { func_name= "strloc"; timer_call,func_name; // Compatibility parameters (including error handling) // compatible.pro: // Parameters for IDL/Yorick compatibility (including error handling) // Array starting index offset, and final index negative offset // (0, 1 in IDL; 1, 0 in Yorick) I0=1 ; I1=0; // Null value ("" in IDL; [] in Yorick) NULL= []; if (is_void(name_list)) return timer_return(func_name, 0); if (is_void(label)) return timer_return(func_name, strloc_err(name_list, "NULL", "", comment=comment)); if (is_integer(label)) { // Integer label index; simply return label index if ((label < 1) || (label > numberof(name_list))) return timer_return(func_name, strloc_err(name_list, strnum(label), "out of range", comment=comment)); return timer_return(func_name, label); } if (typeof(label) != "string") return timer_return(func_name, strloc_err(name_list, typeof(label),"invalid type",comment=comment)); if (label == "") return timer_return(func_name, strloc_err(name_list, label, "", comment=comment)); if (param_set(abbrev)) { // Check for abbreviated match if (param_set(case_fold)) { where_label= where( strtolower(strmid(name_list,0,strlen(label))) == strtolower(label) ); } else { where_label= where( strmid(name_list,0,strlen(label)) == label); } } else { // Check for exact match if (param_set(case_fold)) { where_label= where( strtolower(name_list) == strtolower(label) ); } else { where_label= where( name_list == label ); } } if (numberof(where_label) == 1) return timer_return(func_name, where_label(I0)+I1); if (numberof(where_label) > 1) { // Multiple match; check lengths where_same_len= where(strlen(name_list(where_label)) == strlen(label)); if (numberof(where_same_len) == 1) return timer_return(func_name, where_label(where_same_len(I0))+I1); return timer_return(func_name, strloc_err(name_list, label, "ambiguous", comment=comment)); } return timer_return(func_name, strloc_err(name_list, label, "not found", comment=comment)); } func strloc_err(name_list,label,errmsg,comment=) /* DOCUMENT strloc_err,name_list,label,errmsg,comment= * Error processing for STRLOC. * If COMMENT is null, simply return 0. * If not, display NAME_LIST and ERRMSG, and abort. * SEE ALSO: strloc */ { // Compatibility parameters (including error handling) // compatible.pro: // Parameters for IDL/Yorick compatibility (including error handling) // Array starting index offset, and final index negative offset // (0, 1 in IDL; 1, 0 in Yorick) I0=1 ; I1=0; // Null value ("" in IDL; [] in Yorick) NULL= []; if (is_void(comment)) return 0; if (comment != "") write, comment + "=" ; else write, "Valid choices are:"; nlist= numberof(name_list); for (j=I0; j <= nlist-I1; j++) write, format="%2d%s%s%s\n", j+I1, " or '", name_list(j), "'"; error, "Error: '" + label + "' " + errmsg; } func rangeloc(xcoord,xrange,subscript=,list=,comment=) /* DOCUMENT rangeloc(xcoord,xrange,subscript=0/1,list=0/1,comment=) * Locate nearest value(s) of coordinate range XRANGE on grid XCOORD, and * return the corresponding index value(s), starting from 1. * If XRANGE is a single value, or two equal values, or two unequal values * bracketing a single grid-point, a single index is returned. * Otherwise a two-element array of indices is returned. * If XRANGE is a one-element array, the returned value is also * a one-element array. * If SUBSCRIPT==1, assume that integer coordinate values of XRANGE * represent array subscript values, starting from 1, rather than actual * coordinate values. Subscript values 0, -1, -2, ..., may also be used for * Yorick-style reverse array indices counting from the } of the array. * Of course, floating point values are always assumed to represent actual * coordinate values. * If LIST==1, RANGE contains a list of many coordinate values, * rather than a range of coordinate values. (In this case a monotonically * increasing list of index values is returned.) * XCOORD, XRANGE may also be strings, in which case the COMMENT * parameter is passed to STRLOC when locating XRANGE in XCOORD. * SEE ALSO: strloc, monotonic */ { // Parameters for IDL/Yorick compatibility // Array starting index offset, and final index negative offset // (0, 1 in IDL; 1, 0 in Yorick) I0=1 ; I1=0; // Null value ("" in IDL; [] in Yorick) NULL= ""; if (is_null(xrange) || is_null(xcoord)) return NULL; nr= numberof(xrange); if (param_set(list) && (nr > 1)) { // Handle list option recursively icoord= array(long,nr); for (j=I0; j <= nr-I1; j++) { icoord(j)= rangeloc(xcoord,xrange(j),subscript=subscript,comment=comment); } if (monotonic(icoord) != 1) error, "Coordinate indices not in monotonic order"; return icoord; } subscript_spec= param_set(subscript) && is_integer(xrange); nx= numberof(xcoord); if (subscript_spec) { // Check array subscript specification n0= xrange(I0); if (n0 < 1) n0= nx + n0; if ((n0 < 1) || (n0 > nx)) error, "Subscript out of bounds: " + strnum(xrange(I0)); if (nr == 1) { if (is_scalar(xrange)) return n0 ; else return [ n0 ]; } n1= xrange(I0+1); if (n1 < 1) n1= nx + n1; if ((n1 < 1) || (n1 > nx)) error, "Subscript out of bounds: " + strnum(xrange(I0+1)); // Return index values if (n0 == n1) return [ n0 ]; if (n0 > n1) return [n1, n0] ; else return [n0, n1]; } if (typeof(xcoord) == "string") { // String range; locate first point iloc0= strloc(xcoord, xrange(I0), case_fold=1, abbrev=1, comment=comment); if (iloc0 == 0) error, "Invalid string value - " + xrange(I0); if (nr == 1) { if (is_scalar(xrange)) return iloc0; else return [ iloc0 ]; } // Locate second point iloc1= strloc(xcoord, xrange(I0+1), case_fold=1, abbrev=1, comment=comment); if (iloc1 == 0) error, "Invalid string value - " + xrange(I0+1); if (iloc0 == iloc1) return [ iloc0 ]; if (iloc0 > iloc1) return [iloc1, iloc0]; else return [iloc0, iloc1]; } if (nx == 1) { // Only one coordinate value available if (is_scalar(xrange)) return 1 ; else return [ 1 ]; } // Ascending order flag ascending = (xcoord(I0) < xcoord(nx-I1)); if (nx > 1) { // Check for monotonicity of coordinate values if (ascending) { if (anyof( (xcoord(I0+1:nx-I1) - xcoord(I0:nx-I1-1)) <= 0.) ) error, "Coordinate values not in monotonic order"; } else { if (anyof( (xcoord(I0:nx-I1-1) - xcoord(I0+1:nx-I1)) <= 0.) ) error, "Coordinate values not in monotonic order"; } } // Check for equal values in range nval= nr; if (nr == 2) { if (xrange(I0) == xrange(I0+1)) nval= 1; } if (nval == 1) { // Locate single point x0= xrange(I0); xb= array(double,nx+1); xb(I0)= xcoord(I0); xb(nx+1-I1)= xcoord(nx-I1); xb(I0+1:nx-I1)= 0.5*(xcoord(I0:nx-1-I1) + xcoord(I0+1:nx-I1)); if (ascending) { for (i=I0; i <= nx-I1; i++) if ((x0 > xb(i)) && (x0 <= xb(i+1))) n0= i+I1; if (x0 > xb(nx+1-I1)) n0= nx; if (x0 <= xb(I0)) n0= 1; } else { for (i=I0; i <= nx-I1; i++) if ((x0 <= xb(i)) && (x0 > xb(i+1))) n0= i+I1; if (x0 <= xb(nx+1-I1)) n0= nx; if (x0 > xb(I0)) n0= 1; } if (is_scalar(xrange)) return n0 ; else return [ n0 ]; } else { // Locate range if (ascending) { x0= min(xrange); x1= max(xrange); n0= 1; for (i=I0; i <= nx-I1; i++) if (x0 > xcoord(i)) n0= i+1+I1; if (n0 == nx+1) n0= nx; n1= nx; for (i=nx-I1; i >= I0; i--) if (x1 < xcoord(i)) n1= i-1+I1; if (n1 == 0) n1= 1; } else { x0= max(xrange); x1= min(xrange); n0= 1; for (i=I0; i <= nx-I1; i++) if (x0 < xcoord(i)) n0= i+1+I1; if (n0 == nx+1) n0= nx; n1= nx; for (i=nx-I1; i >= I0; i--) if (x1 > xcoord(i)) n1= i-1+I1; if (n1 == 0) n1= 1; } } if (n0 == n1) return [ n0 ]; if (n0 > n1) return [n1, n0] ; else return [n0, n1]; } func dim_conform( dimsof_a, dimsof_b, trim=, broadcast=) /* DOCUMENT dim_conform(dimsof_a, dimsof_b, trim=trim, broadcast=broadcast) * Returns 1 if the two arrays with dimensions DIMSOF_A and DIMSOF_B are * conformant, i.e., have the same structure. * If TRIM is true, trailing unit-length dimensions in both arrays are ignored. * If BROADCAST==1, unit-length dimensions in array A are allowed to be * broadcast to achieve conformance with array B. * If BROADCAST==2, unit-length dimensions in both arrays are allowed to be * broadcast to achieve conformance. * SEE ALSO: array_eq, dimsof, broadcast, reshape_array */ { func_name= "dim_conform"; timer_call,func_name; // Compatibility parameters I0=1 ; I1=0; if (is_null(dimsof_a) && is_null(dimsof_b)) return timer_return(func_name, 1); if (is_null(dimsof_a) || is_null(dimsof_b)) return timer_return(func_name, 0); ndima= dimsof_a(I0); ndimb= dimsof_b(I0); dima= dimsof_a(I0:I0+ndima); dimb= dimsof_b(I0:I0+ndimb); ndim0= min([ndima, ndimb]); if (param_set(trim) && (ndima != ndimb)) { // Trim trailing unit-length dimensions if (ndima > ndimb) { if (allof(dima(I0+ndim0+1:) == 1)) { // Trim A dima= dima(I0:I0+ndim0); dima(I0)= ndim0; ndima= ndim0; } } else { if (allof(dimb(I0+ndim0+1:) == 1)) { // Trim B dimb= dimb(I0:I0+ndim0); dimb(I0)= ndim0; ndimb= ndim0; } } } if (!param_set(broadcast)) { // No broadcasting if (ndima != ndimb) return timer_return(func_name, 0); return timer_return(func_name, (max(abs(dima - dimb)) == 0)); } else { // Broadcast unit-length dimensions if ((ndima > ndimb) && (broadcast != 2)) return timer_return(func_name, 0); for (i=1; i <= ndim0; i++) { if ( (dima(I0+i) != dimb(I0+i)) && \ (dima(I0+i) != 1) && \ ((dimb(I0+i) != 1) || (broadcast != 2)) ) return timer_return(func_name, 0); } return timer_return(func_name, 1); } } func dim_reshape(dims,insert=,count=,focus=, delete=,trim=,mindim=,broadcast=) /* DOCUMENT dim_reshape(dims,insert=insert,count=count,focus=focus, delete=delete,trim=0/1,mindim=mindim,broadcast= Reshapes a Yorick-style dimension list DIMS. If INSERT >= 0, an extra dimension of length COUNT is inserted to the right of dimension INSERT. (If COUNT is omitted, unit-length is assumed; if there are not INSERT dimensions, extra unit-length dimensions are added so that the result has at least INSERT+1 dimensions.) If FOCUS > 0, a 3-dimensional list is returned, consisting of [3, all dimensions to left, dimension FOCUS, all dimensions to right]. If DELETE > 0, a list with dimension DELETE removed is returned. (If there are not DELETE dimensions, nothing is done; only unit-length dimensions may be deleted.) If TRIM == 1, non-unit length trailing dimensions are removed. If MINDIM > 0, at least MINDIM dimensions are guaranteed on the list, through padding with unit-length dimensions. (The above options may be combined, in the above order) If BROADCAST is set to Yorick-style dimension list, then DIMS is re-shaped to achieve broadcast conformance between both dimension lists. SEE ALSO: dimsof, dim_conform, broadcast */ { func_name= "dim_reshape"; timer_call,func_name; // Compatibility parameters I0=1; I1=0 dim0= dims; ndim0= dim0(I0); if (!is_null(broadcast)) { dim1= broadcast; ndim1= dim1(I0); ndim2= max([ndim0, ndim1]); dim2= array(1, ndim2+1); dim2(I0)= ndim2; for (i=I0+1; i <= I0+ndim2; i++) { if ((i <= I0+ndim0) && (i <= I0+ndim1)) { if ((dim0(i) != dim1(i)) && (dim0(i) != 1) && (dim1(i) != 1)) error, "Dimensions are not broadcast conformable"; dim2(i)= max([dim0(i), dim1(i)]); } else { if (i <= ndim0) dim2(i)= dim0(i); if (i <= ndim1) dim2(i)= dim1(i); } } return timer_return(func_name, dim2); } if (!is_null(insert)) { newdims= array(long(1), max([insert+2,ndim0+2])); newdims(I0)= numberof(newdims)-1; ndim1= min([ndim0+1, insert+1]); if (ndim1 > 1) newdims(I0+1:I0+ndim1-1)= dim0(I0+1:I0+ndim1-1); if (insert < ndim0) newdims(I0+insert+2:I0+ndim0+1)= dim0(I0+insert+1:I0+ndim0); if (param_set(count)) newdims(I0+insert+1)= count; dim0= newdims; ndim0= newdims(I0); } if (param_set(focus)) { // Create dimension lengths (extended with unit dimensions) len= array(long(1), 2+max([ndim0, focus]) ); if (ndim0 > 0) len(I0+1:I0+ndim0)= dim0(I0+1:I0+ndim0); dim0= [3, prod(len(I0:I0+focus-1)), len(I0+focus), prod(len(I0+focus+1:)) ]; ndim0= 3; } if (param_set(delete)) { if (delete <= ndim0) { if (dim0(I0+delete) != 1) error, "Cannot delete non-unit-length dimensions"; newdims= ndim0-1; if (delete > 1) grow, newdims, dim0(I0+1:I0+delete-1); if (delete < ndim0) grow, newdims, dim0(I0+delete+1:); dim0= newdims; ndim= newdims(I0); }; } if (param_set(trim)) { // Trim trailing unit-length dimensions lastdim= 0; for (i=I0+1; i<=I0+ndim0; i++) { if (dim0(i) > 1) lastdim= i-I0; } newdims= array(long(lastdim), lastdim+1); if (lastdim > 0) newdims(I0+1:I0+lastdim)= dim0(I0+1:I0+lastdim); dim0= newdims; ndim0= newdims(I0); } if (param_set(mindim)) { if (ndim0 < mindim) { // Extend to minimum number of dimensions newdims= array(long(1), mindim+1); newdims(I0)= mindim; if (ndim0 > 0) newdims(I0+1:I0+ndim0)= dim0(I0+1:I0+ndim0); dim0= newdims; ndim0= newdims(I0); } } return timer_return(func_name, dim0); } func array_eq( array_a, array_b, trim=, epsilon=) /* DOCUMENT array_eq(array_a, array_b, trim=, epsilon=) * is true if the two arrays ARRAY_A and ARRAY_B are identical. * If TRIM is true, trailing unit-length dimensions in both arrays are ignored. * If EPSILON > 0, it represents the allowed fractional difference * (expressed a fraction of the maximum absolute value among both arrays). * SEE ALSO: dim_conform */ { func_name= "array_eq"; timer_call,func_name; if (is_null(array_a) && is_null(array_b)) return timer_return(func_name, 1); if (is_null(array_a) || is_null(array_b)) return timer_return(func_name, 0); if (!dim_conform(dimsof(array_a), dimsof(array_b), trim=trim)) return timer_return(func_name, 0); if (param_set(epsilon) && is_number(array_a)) { // Specified tolerance maxabs= abs(epsilon)*max([max(abs(array_a)), max(abs(array_b))]); return timer_return(func_name, allof( abs(array_a(*)-array_b(*)) <= maxabs )); } else { // Zero tolerance return timer_return(func_name,allof(array_a(*) == array_b(*))); } } func alleq( arr, epsilon=) /* DOCUMENT alleq(arr, epsilon=) * is true if all elements of array ARR are equal to each other. * If EPSILON > 0, it represents the allowed fractional difference * (expressed a fraction of the maximum absolute value in the array). * SEE ALSO: allof, noneof, anyof, monotonic */ { func_name= "alleq"; timer_call,func_name; if (param_set(epsilon) && is_number(arr)) { // Specified tolerance maxabs= abs(epsilon)*max(abs(arr)); return timer_return(func_name, allof( abs(arr(*)-arr(1)) <= maxabs )); } else { // Zero tolerance return timer_return(func_name,allof( arr(*) == arr(1) )); } } func monotonic(arr) /* DOCUMENT monotonic(arr) * Return 1 if ARR is in strictly ascending order, * -1 if it is in strictly descending order, * or 0 if it is non-monotonic or has equal values, or has a single value. * SEE ALSO: rangeloc, alleq */ { // Compatibility parameters (including error handling) // compatible.pro: // Parameters for IDL/Yorick compatibility (including error handling) // Array starting index offset, and final index negative offset // (0, 1 in IDL; 1, 0 in Yorick) I0=1 ; I1=0; // Null value ("" in IDL; [] in Yorick) NULL= []; // Number of coordinate values nx= numberof(arr); if (nx == 1) return 0; // Check for monotonicity of coordinate values ret_code= 0; if (arr(I0) < arr(I0+1)) { // Check for ascending order if (allof( (arr(I0+1:nx-I1) - arr(I0:nx-I1-1)) > 0) ) ret_code= 1; } else { // Check for descending order if (allof( (arr(I0:nx-I1-1) - arr(I0+1:nx-I1)) > 0) ) ret_code= -1; } return ret_code; } func reshape_array( &arr, dims) /* DOCUMENT reshape_array, arr, dims Reshapes array ARR to have dimensions DIMS SEE ALSO: array, dim_reshape, rangeop, broadcast, grow */ { func_name= "reshape_array"; timer_call,func_name; if ( dim_conform(dimsof(arr),dims) ) return timer_return(func_name); if (dims(1) == 0) { // Scalar reshape if (numberof(arr) != 1) error, "Incorrect size for reshape"; } else { // Array reshape if (numberof(arr) != prod(dims(2:))) error, "Incorrect size for reshape"; } // Reshape array temarr= array(structof(arr), dims) temarr(*)= arr(*) arr= temarr ///reshape, arr, structof(arr), dims; return timer_return(func_name); } func xrotate(xperiod,rot_angle, &xcoord,&rot_index,&xrange) //YORICKoutput: /* DOCUMENT xrotate,xperiod,rot_angle,xcoord,rot_index,xrange * Carries out rotation of periodic coordinate dimension * Input parameters * Input: xperiod -- coordinate period (e.g., 360 for degrees) * rot_angle -- rotation angle (in the same units as the coordinate) * Input/output: * xcoord -- X-coordinate vector to be rotated * Output: * rot_index -- rotation index * (i.e, subscript of new first element, starting from 0) * Input/output: * xrange -- X-coordnate range ([xmin,xmax]) to be rotated (optional) * * SEE ALSO: rangeop */ { func_name= "xrotate"; timer_call,func_name; // Parameters for IDL/Yorick compatibility // Array starting index offset, and final index negative offset // (0, 1 in IDL; 1, 0 in Yorick) I0=1 ; I1=0; // Null value ("" in IDL; [] in Yorick) NULL= []; // Ensure that 0 <= rot_angle < xperiod angle= rot_angle; while (angle < 0) angle= angle + xperiod; while (angle > xperiod) angle= angle - xperiod; // Coordinate count nx= numberof(xcoord); // Reversed coordinate flag rev_flag= xcoord(nx-I1) < xcoord(I0); if (rev_flag) xcoord= rangeop(xcoord,"rev",1); // New origin index org_index= rangeloc( xcoord, xcoord(I0)+angle) - 1; // Initialize rotated X coordinate vector new_xcoord= xcoord(I0+org_index:); shift_flag= 0; if (org_index > 0) { // Append wraparound segment of X coordinate grow, new_xcoord, xcoord(I0:org_index-I1)+xperiod; // Coordinate shift flag shift_flag= (new_xcoord(I0) >= 0.499*xperiod); if (shift_flag) new_xcoord= new_xcoord - xperiod; } // Rotation index rot_index= org_index; if (rev_flag) { // Un-reverse coordinate new_xcoord= rangeop(new_xcoord,"rev",1); rot_index= (nx - rot_index) % nx; } // Copy rotated X-coordinate xcoord= new_xcoord; if (!is_null(xrange)) { // Rotate X range xrange= xrange + angle; if (shift_flag) xrange= xrange - xperiod; } return timer_return(func_name); } func broadcast( arr, dims) /* DOCUMENT broadcast(arr,dims) * Carries out Yorick-style broadcasting of array ARR to conform to * dimension list DIMS. * SEE ALSO: rangeop, grow, reshape_array, outer_prod */ { func_name= "broadcast"; timer_call,func_name; // Parameters for IDL/Yorick compatibility // Array starting index offset, and final index negative offset I0=1 ; I1=0; // Null value NULL= []; if ( dim_conform(dimsof(arr),dims) ) return timer_return(func_name,arr); // Number of dimensions dima= dimsof(arr); ndima= dima(I0); ndimb= dims(I0); ndim0= min([ndima, ndimb]); ndim1= max([ndima, ndimb]); temarr= arr; for (i=1; i <= ndim0; i++) { if (dima(I0+i) != dims(I0+i)) { // Dimension lengths differ if (dima(I0+i) == 1) { // Broadcast dimension temarr= rangeop( temarr, "bcst", i, count=dims(I0+i), temporary=1 ); } else { error, "Dimension length mismatch"; } } } if (ndima != ndimb) { if (ndima > ndimb) { if (anyof(dima(I0+ndim0+1:) > 1)) error, "Cannot truncate non-unit length dimensions"; // Truncate array by reshaping reshape_array, temarr, dims; } else { // Broadcast extra dimensions temarr= rangeop( temarr, "-", count=dims(I0+ndim0+1:), temporary=1 ); } } return timer_return(func_name, temarr); } func rangeop(arr, range, mdim, count=, missing_value=, pad=, keep=, temporary=) /* DOCUMENT rangeop(arr, range, mdim, count=, missing_value=, * pad=, keep=0/1, temporary=0/1 ) * Carries out Yorick-style range operation RANGE * on dimension MDIM of array ARR. * MDIM defaults to the last dimension, if omitted. * RANGE=index,[index1,index2],"-"/"bcst"/"swap"/"shift"/"rot"/"rev"/"sum"/"min"/"max") * (For "sum"/"min"/"max" operations, if MISSING_VALUE is specified, * missing values are ignored during the operation.) * RANGE == "bcst" broadcasts an unit-length dimension to have length COUNT. * RANGE == "swap" swaps two halves of a dimension * RANGE == "shift" shifts dimension by COUNT with no wraparound * (PAD=padding_value may be specified for padded shifts) * RANGE == "rot" rotates dimension by COUNT with wraparound * For "shift"/"rot", COUNT > 0 moves elements rightward (forward) by COUNT, * and COUNT < 0 moves elements leftward by COUNT; * RANGE == "rev" reverses a dimension * RANGE == "-" inserts new dimension(s) to the right of MDIM, where MDIM * may range from 0 to the total number of dimensions * and COUNT is the length of the dimension(s) to be inserted * (COUNT may be an array) * If KEEP==1, keep all dimensions during averaging operations * i.e., by reducing averaged dimensions to unit-length * If TEMPORARY==1, input array is assumed to be a temporary array and may * be overwritten * SEE ALSO: arrayop, xrotate, broadcast, grow, reshape_array */ { func_name= "rangeop"; timer_call,func_name; // Parameters for IDL/Yorick compatibility // Array starting index offset, and final index negative offset I0=1 ; I1=0; // Null value NULL= []; if (is_null(arr)) error, "Error - null array"; // Number of dimensions adims= dimsof(arr); ndim0= adims(I0); //IDLbegin: // Pseudo-index flag //:pseudoind= 0; //:if (typeof(range) == "string") pseudoind= (range == "-"); //:if ((ndim0 == 0) && (!pseudoind)) { // Pretend scalars are 1-D arrays (for non-pseudo-index operations) //: adims= [1, 1]; //: ndim0= 1; //:} //IDLend: if (is_null(mdim)) mdim= ndim0; if (typeof(range) != "string") { // Index range if ((!is_integer(range)) || (numberof(range) > 2)) error, "Invalid range operation"; if (mdim < 1) error, "Invalid dimension number"; tem_dims= dim_reshape(adims,focus=mdim); leftdim= tem_dims(I0+1); middim= tem_dims(I0+2); rightdim= tem_dims(I0+3); if (middim == 1) { // Unit length dimension; accept any range range2= [I0, I0]; } else { if (numberof(range) == 1) { // Single subscript if ((range(I0) < I0) || (range(I0) > middim-I1)) error, "Array subscript out of range"; range2= [range(I0), range(I0)]; } else { // Subscript range if ( (range(I0) < I0) || (range(I0+1) > middim-I1) || (range(I0) > range(I0+1)) ) error, "Array subscripts out of range"; range2= range; } } // New count in selected dimension new_dims= dim_reshape( adims, mindim=mdim ); new_dims(I0+mdim)= range2(I0+1) - range2(I0) + 1; // Reshape array to have 3 dimensions reshape_array, arr, tem_dims; // Create 3-D output array of same type new_arr= array(arr(I0), dim_reshape(new_dims,focus=mdim) ); // Select subscript range new_arr(,,)= arr(,range2(I0):range2(I0+1),); } else { // String range if (range == "bcst") { // Broadcast dimension if (mdim < 1) error, "Invalid dimension number"; tem_dims= dim_reshape(adims,focus=mdim); leftdim= tem_dims(I0+1); middim= tem_dims(I0+2); rightdim= tem_dims(I0+3); if (middim != 1) error, "Cannot broadcast non-unit length dimension"; if (is_null(count)) error, "COUNT must be specified for broadcasting dimension"; if (!is_scalar(count)) error, "COUNT parameter must be a scalar"; // New dimensions; copy broadcast dimension count if (mdim <= ndim0) { new_dims= adims; } else { new_dims= array(1,mdim+1); new_dims(I0)= mdim; new_dims(I0+1:I0+ndim0)= adims(I0+1:I0+ndim0); } new_dims(I0+mdim)= count; // Reshape array to have 2 dimensions reshape_array, arr, [2, leftdim, rightdim]; // Broadcast dimension new_arr= arr(,-:1:count,) } else if ((range == "swap") || (range == "shift") || (range == "rot") || (range == "rev")) { // Swap/shift/rotate/reverse dimension if (mdim < 1) error, "Invalid dimension number"; tem_dims= dim_reshape(adims,focus=mdim); leftdim= tem_dims(I0+1); middim= tem_dims(I0+2); rightdim= tem_dims(I0+3); // New dimensions new_dims= adims; // Reshape array to have 3 dimensions reshape_array, arr, tem_dims; // Create 3-D output array of same type new_arr= array(arr(I0), dim_reshape(new_dims,focus=mdim) ); if (range == "swap") { // Swap two halves of dimension ntot= adims(I0+mdim); nhalf= ntot/2; if ((ntot % 2) == 1) new_arr(,I0+nhalf,)= arr(,I0+nhalf,); if (ntot > 1) { new_arr(,I0:I0+nhalf-1,)= arr(,ntot-I1-nhalf+1:,); new_arr(,ntot-I1-nhalf+1:,)= arr(,I0:I0+nhalf-1,); } } else if (range == "shift") { // Shift dimension ntot= adims(I0+mdim); if (!param_set(count)) { // No shift performed new_arr(*)= arr(*); } else { if (abs(count) >= ntot) error, "Shift count out of bounds"; if (count > 0) { // Right (forward) shift new_arr(,I0+count:,)= arr(,I0:ntot-I1-count,); if (is_null(pad)) { new_arr(,I0:I0+count-1,)= (arr(,I0,))(,-:1:count,); } else { new_arr(,I0:I0+count-1,)= pad; } } else { // Left (backward) shift bcount= -count; new_arr(,I0:ntot-I1-bcount,)= arr(,I0+bcount:,); if (is_null(pad)) { new_arr(,ntot-I1-bcount+1:,)= (arr(,ntot-I1,))(,-:1:bcount,); } else { new_arr(,ntot-I1-bcount+1:,)= pad; } } } } else if (range == "rot") { // Rotate dimension ntot= adims(I0+mdim); if (!param_set(count)) rcount= 0 ; else rcount= count % ntot; if (rcount == 0) { // No rotation performed new_arr(*)= arr(*); } else { if (rcount > 0) { // Right (forward) rotation new_arr(,I0+rcount:,)= arr(,I0:ntot-I1-rcount,); new_arr(,I0:I0+rcount-1,)= arr(,ntot-I1-rcount+1:,); } else { // Left (backward) rotation bcount= -rcount; new_arr(,I0:ntot-I1-bcount,)= arr(,I0+bcount:,); new_arr(,ntot-I1-bcount+1:,)= arr(,I0:I0+bcount-1,); } } } else if (range == "rev") { // Reverse dimension new_arr(,,)= arr(,::-1,); //IDL2YORICK: new_arr= reverse(arr,2) } else error, "Invalid case option"; } else if (range == "-") { // Pseudo-index if (mdim < 0) error, "Invalid dimension number"; // New count in each dimension if (is_null(count)) { // Default: unit-length dimension countlist= 1; } else { // Copy dimension lengths to be inserted countlist= count; } ncount= numberof(countlist); new_dims= array(long(1), max([ndim0,mdim]) + ncount + 1 ); new_dims(I0)= numberof(new_dims)-1; new_dims(I0+mdim+1:I0+mdim+ncount)= countlist(*); if (mdim > 0) new_dims(I0+1:I0+mdim)= adims(I0+1:I0+mdim); if (mdim < ndim0) new_dims(I0+mdim+ncount+1:)= adims(I0+mdim+1:); // Focus array tem_dims= dim_reshape(adims,insert=mdim,focus=mdim+1); leftdim= tem_dims(I0+1); rightdim= tem_dims(I0+3); // Reshape array to have 2 dimensions reshape_array, arr, [2, leftdim, rightdim]; // Create 3-D output array of same type new_arr= array(arr(I0), [3, leftdim, prod(countlist), rightdim]); // Pseudo-index: insert dimension new_arr(,,)= arr(,-,); } else { if (mdim < 1) error, "Invalid dimension number"; // If scalar, simply return scalar if (ndim0 == 0) return timer_return(func_name, arr(I0)); tem_dims= dim_reshape(adims,focus=mdim); leftdim= tem_dims(I0+1); middim= tem_dims(I0+2); rightdim= tem_dims(I0+3); // New dimensions (make averaged dimension of unit length) new_dims= dim_reshape( adims, mindim=mdim ); new_dims(I0+mdim)= 1; // Reshape input array to have 3 dimensions reshape_array, arr, tem_dims; // Create 3-D output array of same type new_arr= array(arr(I0), dim_reshape(new_dims,focus=mdim) ); if (!is_null(missing_value)) { // Locate missing values in input array where_missing= where(arr == missing_value); if (is_where(where_missing)) { if (range == "sum") { arr(where_missing)= 0.; } else if (range == "min") { // Substitute missing values with maximum array value maxarr= max(arr); arr(where_missing)= maxarr; } else if (range == "max") { // Substitute missing values with minimum array value minarr= min(arr); arr(where_missing)= minarr; } } } else { where_missing= NULL; } // Carry out requested range operation if (range == "sum") { // Find sum along dimension new_arr(,I0,)= arr(,sum,); } else if (range == "min") { // Find minimum along dimension new_arr(,I0,)= arr(,min,); } else if (range == "max") { // Find maximum along dimension new_arr(,I0,)= arr(,max,); } else { error, "Invalid range operation '" + range + "'"; } // Put missing values back if input array is not temporary if (is_where(where_missing) && !param_set(temporary)) arr(where_missing)= missing_value; if (!param_set(keep)) { // Remove "averaged" dimension new_dims= dim_reshape(new_dims, delete=mdim); } } } // Reshape input array to original shape, if not temporary if (!param_set(temporary)) reshape_array, arr, adims; // Reshape output array to final shape reshape_array, new_arr, new_dims; return timer_return(func_name, new_arr); } func arrayop(arr,range,mdim,missing_value=,list=) /* DOCUMENT arrayop(arr, range, mdim, missing_value=, list=0/1) * Carries out Yorick-style range operation RANGE * on dimension MDIM of array ARR. * (Efficient version of RANGEOP for up to five dimensions) * MDIM defaults to the last dimension, if omitted. * NOTE: A "reduced" dimension is not eliminated, but is retained as an unit * length dimension (unlike in RANGEOP). * RANGE=index,[ index1,index2], "sum"/"min"/"max" ) * (For "sum"/"min"/"max" operations, if MISSING_VALUE is specified, * missing values are ignored during the operation.) * If LIST==1, RANGE contains a list of array indices for dimension MDIM, * rather than a range of indices. * SEE ALSO: rangeop */ { func_name= "arrayop"; timer_call,func_name; // Compatibility parameters (including error handling) // compatible.pro: // Parameters for IDL/Yorick compatibility (including error handling) // Array starting index offset, and final index negative offset // (0, 1 in IDL; 1, 0 in Yorick) I0=1 ; I1=0; // Null value ("" in IDL; [] in Yorick) NULL= []; // Number of dimensions adims= dimsof(arr); ndim0= adims(I0); if (is_null(mdim)) mdim= max([1, ndim0]); if (ndim0 > 5) error, "Cannot handle more than 5 dimensions"; if ((mdim < 1) || (mdim > 5)) error, "Invalid dimension number " + strnum(mdim); if (typeof(range) == "string") { // String range if (!is_null(missing_value)) return timer_return(func_name, rangeop( arr, range, mdim, missing_value=missing_value, keep=1 )); if (ndim0 == 0) return timer_return(func_name, arr); new_count= (dim_reshape(adims, mindim=mdim))(I0+1:); new_count(mdim-I1)= 1; if (range == "sum") { // Find sum along dimension //IDLbegin: //: return timer_return(func_name, reform( sum(arr,mdim), new_count )); //IDLend: //YORICKbegin: if (mdim == 1) { return timer_return(func_name,arr(sum,..)(-,..)); } else if (mdim == 2) { return timer_return(func_name,arr(,sum,..)(,-,..)); } else if (mdim == 3) { return timer_return(func_name,arr(,,sum,..)(,,-,..)); } else if (mdim == 4) { return timer_return(func_name,arr(,,,sum,..)(,,,-,..)); } else if (mdim == 5) { return timer_return(func_name,arr(,,,,sum,..)(,,,,-,..)); } //YORICKend: } else if (range == "min") { // Find minimum along dimension //IDLbegin: //: return timer_return(func_name, timer_return(func_name, rangeop( arr,range,mdim,missing_value=missing_value,keep=1 ))); //IDLend: //YORICKbegin: if (mdim == 1) { return timer_return(func_name,arr(min,..)(-,..)); } else if (mdim == 2) { return timer_return(func_name,arr(,min,..)(,-,..)); } else if (mdim == 3) { return timer_return(func_name,arr(,,min,..)(,,-,..)); } else if (mdim == 4) { return timer_return(func_name,arr(,,,min,..)(,,,-,..)); } else if (mdim == 5) { return timer_return(func_name,arr(,,,,min,..)(,,,,-,..)); } //YORICKend: } else if (range == "max") { // Find maximum along dimension //IDLbegin: //: return timer_return(func_name, timer_return(func_name, rangeop( arr,range,mdim,missing_value=missing_value,keep=1 ))); //IDLend: //YORICKbegin: if (mdim == 1) { return timer_return(func_name,arr(max,..)(-,..)); } else if (mdim == 2) { return timer_return(func_name,arr(,max,..)(,-,..)); } else if (mdim == 3) { return timer_return(func_name,arr(,,max,..)(,,-,..)); } else if (mdim == 4) { return timer_return(func_name,arr(,,,max,..)(,,,-,..)); } else if (mdim == 5) { return timer_return(func_name,arr(,,,,max,..)(,,,,-,..)); } //YORICKend: } else { error, "Invalid range operation - " + range; } } else { // Index range/list if (param_set(list) && (numberof(range) > 1)) { // Index list if (ndim0 == 0) error, "Invalid index list for scalar"; if (mdim == 1) { return timer_return(func_name, arr(range,,,,)); } else if (mdim == 2) { return timer_return(func_name, arr(,range,,,)); } else if (mdim == 3) { return timer_return(func_name, arr(,,range,,)); } else if (mdim == 4) { return timer_return(func_name, arr(,,,range,)); } else if (mdim == 5) { return timer_return(func_name, arr(,,,,range)); } } else { // Index range if (numberof(range) == 1) { // Single index irange= [range(I0), range(I0)]; } else { // Index range irange= range; } if (ndim0 == 0) { if ( (irange(I0) != I0) || (irange(I0) != irange(I0+1)) || \ (numberof(range) > 2) ) error, "Invalid index for scalar"; return timer_return(func_name, arr); } if (mdim == 1) { return timer_return(func_name, arr(irange(I0):irange(I0+1),,,,)); } else if (mdim == 2) { return timer_return(func_name, arr(,irange(I0):irange(I0+1),,,)); } else if (mdim == 3) { return timer_return(func_name, arr(,,irange(I0):irange(I0+1),,)); } else if (mdim == 4) { return timer_return(func_name, arr(,,,irange(I0):irange(I0+1),)); } else if (mdim == 5) { return timer_return(func_name, arr(,,,,irange(I0):irange(I0+1))); } } } } func lonlabel( axis, index, value) /* DOCUMENT lonlabel(axis,index,value) * Returns longitude coordinate label E or W, * given the latitude VALUE in degrees North * (AXIS and INDEX parameters are ignored, and are present for * compatibility with IDL tick formats) * SEE ALSO: hlegend, lonlabel, datelabel */ { format= "%7.2f" //IDL2YORICK: format= "(F7.2)" ; lon = value; if (lon > 180) lon= lon - 360; if (lon > 0) return strnum(lon,format=format)+"E"; if (lon < 0) return strnum(-lon,format=format)+"W"; if (lon == 0) return "0"; if (lon == 180) return "180"; } func latlabel( axis, index, value) /* DOCUMENT latlabel(axis,index,value) * Returns latitude coordinate label N or S, * given the latitude VALUE in degrees North * (AXIS and INDEX parameters are ignored, and are present for * compatibility with IDL tick formats) * SEE ALSO: hlegend, lonlabel, datelabel */ { format= "%6.2f" //IDL2YORICK: format= "(F6.2)" ; if (value > 0) return strnum(value,format=format)+"N"; if (value < 0) return strnum(-value,format=format)+"S"; return "EQ"; } func datelabel( axis, index, value) /* DOCUMENT datelabel(axis,index,value) * Returns date coordinate labels yyyy/mm/dd * given the date VALUE yyyymmdd.fff * (assuming 0000.000 value for the first day of the first month) * (AXIS and INDEX parameters are ignored, and are present for * compatibility with IDL tick formats) * SEE ALSO: hlegend, lonlabel, latlabel */ { yr= long(value / double(10000.)); rem= value - yr*double(10000.); mm= long(rem / double(100.) ); rem= rem - mm*double(100.); dd= long( rem ); rem= rem - dd; // Correct for zero offset mm= mm+1; dd= dd+1; // Date string datestr= strnum(yr) + "/" + strnum(mm) + "/" + strnum(dd); if (rem >= double(0.001)) { // Include fractional part of day (rounded to three decimal places) ifrac= long(rem/double(0.001) + 0.5); if (ifrac > 0) datestr= datestr + "." +strmid(strnum(1000+ifrac),1,3); } return datestr; } func param_set(var) /* DOCUMENT param_set(var) is true if VAR is not NULL, is not the scalar number zero, is not the null string, and is not a range. SEE ALSO: is_null, is_number, is_scalar */ { vartype= typeof(var); if ((vartype == "void") || (vartype == "range")) return 0; if (is_array(var) && !dimsof(var)(1)) { // Scalar if (vartype == "pointer") return (var != pointer(0) ); if (vartype == "string") return (var != ""); if (vartype == "struct_instance") return 1; return (var != 0); } return 1; } func imaginary(var) /* DOCUMENT imaginary(var) returns imaginary part of complex variable VAR SEE ALSO: */ { if (typeof(var) != "complex") error, "Cannot obtain imaginary part of type "+typeof(var); return var.im; } func is_null(var) /* DOCUMENT is_null(var) * is true if VAR is void, or a null pointer * SEE_ALSO: is_integer, is_number, is_array, typeof */ { vartype= typeof(var); return (vartype == "void") || ( (vartype == "pointer") && !dimsof(var)(1) && (var == pointer(0)) ); } func is_integer(var) /* DOCUMENT is_integer(var) is true if VAR is of an integer type: char/short/int/long SEE ALSO: is_null, is_number, is_scalar */ { return (anyof(strmatch(["char", "short", "int", "long"],typeof(var)))); } func is_number(var) /* DOCUMENT is_number(var) is true if VAR is of an integer or floating-point type: char/short/int/long/float/double SEE ALSO: is_null, is_integer, is_scalar */ { return (anyof(strmatch(["char", "short", "int", "long", "float", "double"], typeof(var)))); } func is_where(where_output) /* DOCUMENT is_where(where_output) is true if WHERE_OUTPUT is non-null range (in Yorick) or is a numeric array of one or more dimensions (in IDL). SEE ALSO: is_null, is_scalar */ { return (numberof(where_output) > 0); } func prod(arr) /* DOCUMENT prod(arr) Returns the product of all elements of an array. SEE_ALSO: sum */ { val= long(1); for (i=1; i <= numberof(arr); i++) val= val * arr(i); return val; } func outer_prod( vec1, vec2) /* DOCUMENT outer_prod(vec1, vec2) * Computes outer product of two vectors VEC1 and VEC2, * returning a matrix with dimensions numberof(VEC1)xnumberof(VEC2). * (i.e., VEC1 # transpose(VEC2) in IDL; VEC1 * VEC2(-,) in Yorick) * SEE ALSO: broadcast */ { func_name= "outer_prod"; timer_call,func_name; return timer_return(func_name,vec1 * vec2(-,)); } func randomu1( &seed, //YORICKoutput: dims=, time_seed=) /* DOCUMENT randumu1(seed, dims=, time_seed=0/1) * Returns a uniformly-distributed pseudo-random number in the interval * from 0.0 to 1.0 in Yorick (always < 1.0 in IDL). * The returned value type is double (promoted from float in IDL). * If DIMS is set to dimension list, an array of random numbers with dimension * DIMS is returned. * In IDL, SEED is an integer. * In Yorick, SEED is a double value in the open interval 0 < SEED < 1. and * If TIME_SEED==1, a new seed is generated using the system clock, and its * value is returned as SEED. * Is SEED is omitted, a continuing sequence of random numbers is generated * using the previous seed (or the default seed at the very first time). * If SEED is specified, a new sequence of random numbers is initialized using * SEED. * SEE ALSO: randomu, randomn */ { func_name= "randomu1"; timer_call,func_name; if (param_set(time_seed)) { // Create new seed using current date and time fmt="%3s %3s %2d %2d:%2d:%2d %4s"; daystr=monstr=yearstr= ""; dd=hh=mm=ss= 0; if (sread(timestamp(),format=fmt,daystr,monstr,dd,hh,mm,ss,yearstr) != 7) error, "Error in reading TIMESTAMP output"; digit12= long((100.*ss)/62.); digit34= long((100.*mm)/61.); digit7= long((10.*hh)/24.); digit8= long((10.*(dd-1))/31.); elapsed= array(0., 3); timer, elapsed; seed= digit12/1.0e2 + (digit34 + sum(abs(elapsed))%1.0)/1.0e4 + digit7/1.0e7 + digit8/1.0e8; // print, daystr,monstr,dd,hh,mm,ss,yearstr; // print, digit12, digit34, sum(abs(elapsed))%1.0, digit7, digit8; // print, seed; } if (!is_void(seed)) { // Re-seed if (is_integer(seed) || (seed <= 0.0) || (seed >= 1.0)) error, "SEED should be a real number > 0, and < 1"; random_seed, seed; } return timer_return(func_name,random(dims)); } func low_pass_filt(m, delt, freq) /* DOCUMENT low_pass_filt(m, delt, freq) * Returns weights of a low-pass finite impulse response filter, * of order 2*M+1, given the sampling interval DELT (in sec), and the * the cut-off frequency FREQ (in Hz) * To be safe, one might choose M >= DELT/FREQ. * (Taken from "Applied Time Series Analysis: Vol. I", * by R.K. Otnes and L. Enochson, p. 169) * SEE ALSO: */ { // Compatibility parameters (including error handling) // compatible.pro: // Parameters for IDL/Yorick compatibility (including error handling) // Array starting index offset, and final index negative offset // (0, 1 in IDL; 1, 0 in Yorick) I0=1 ; I1=0; // Null value ("" in IDL; [] in Yorick) NULL= []; // Magic numbers d0 = 0.35577019 //IDL2YORICK: d0 = 0.35577019D0 ; d3 = [ 0.2436983, 0.07211497, 0.00630165 ] //IDL2YORICK: d3 = [ 0.2436983D0, 0.07211497D0, 0.00630165D0 ] ; d = d0; grow, d, 2.0*d3; // First generate plain boxcar weights wt = array(double,m+1); fact = double(2.0) * freq * delt; wt(I0) = fact; fact = fact * pi(); for (j=1; j <= m; j++) wt(I0+j) = sin(fact*j)/(pi()*j); // Trapezoidal weighting at } wt(m) = wt(m) / double(2.0); // Now apply the Potter P310 window sumg = wt(I0); for (j=1; j <= m; j++) { fact = (pi()*j)/double(m); win = sum( d * cos(fact*(indgen(4)-I0)) ); wt(I0+j) = wt(I0+j) * win; sumg = sumg + double(2.0)*wt(I0+j); } // Normalized weights wt= wt/sumg; // Reflect weights to create symmetric filter full_wt= rangeop(wt(I0+1:),"rev"); grow, full_wt, wt; return full_wt; } func ref(var,noscalar=) /* DOCUMENT ref(var,noscalar=) Returns reference to ("address of") variable VAR, or null, if IS_VOID(VAR) is true. If NOSCALAR==1, a null value is returned if IS_SCALAR(VAR) is true. (Extension of operator "&" Yorick-style; no-op in IDL) SEE ALSO: is_scalar */ { if (param_set(noscalar)) { if (is_scalar(var) && !(typeof(var) == "struct_instance")) return pointer([]); } if (is_null(var)) return pointer([]); return &var; } func deref(ptr) /* DOCUMENT deref(ptr) Deference pointer PTR, but return null value if PTR is itself null. (Error-free version of Yorick-style "*" operator; no-op in IDL) SEE ALSO: ref, * */ { // Parameters for IDL/Yorick compatibility // Array starting index offset, and final index negative offset // (0, 1 in IDL; 1, 0 in Yorick) I0=1 ; I1=0; // Null value ("" in IDL; [] in Yorick) NULL= []; if (is_null(ptr)) return NULL; return *ptr; } func typeconv(typestr,var) /* DOCUMENT typeconv(typestr,var) * Convert variable VAR to type TYPESTR, and return the converted variable * TYPESTR may be one of the following: * "char", "int", "long", "float", "double", "complex", "string" * SEE ALSO: typeof, array */ { func_name= "typeconv"; timer_call,func_name; // Parameters for IDL/Yorick compatibility // Array starting index offset, and final index negative offset I0=1 ; I1=0; // Null value NULL= []; str= strtolower(typestr); if (str == "char") { return timer_return(func_name, char(var)); } else if (str == "int") { return timer_return(func_name, int(var)); } else if (str == "long") { return timer_return(func_name, long(var)); } else if (str == "float") { return timer_return(func_name, float(var)); } else if (str == "double") { return timer_return(func_name, double(var)); } else if (str == "complex") { return timer_return(func_name, complex(var)); } else if (str == "string") { return timer_return(func_name, string(var)); } else { error, "Invalid type '" + str + "'"; } } func raw_dump( &arr, &fname_or_handle, //YORICKoutput: byte_offset=, byte_skip=, rec_cluster=, input=, close_file=) /* DOCUMENT raw_dump, arr, fname_or_handle, byte_offset=, byte_skip=, * rec_cluster=, input=0/1, close_file=0/1 * Write raw data from array ARR to file FNAME_OR_HANDLE, skipping * BYTE_OFFSET (default 0) bytes from the beginning of the file. * If FNAME_OR_HANDLE contains a file name, that file opened, and the * corresponding file handle is returned on output. * If BYTE_SKIP > 0, ARR is split into records, using the last dimension * as the record dimension, with BYTE_SKIP denoting the separation between * clusters of REC_CLUSTER (default 1) records each. * (If ARR is a scalar, it is assumed to represent a single record.) * * If INPUT==1, read data into array ARR from file handle FHANDLE. * * If End-of-File is reached, the ARR is set to null value on return. * * If CLOSE_FILE==1, the file is closed prior to returning. * * SEE ALSO: */ { // Compatibility parameters (including error handling) // compatible.pro: // Parameters for IDL/Yorick compatibility (including error handling) // Array starting index offset, and final index negative offset // (0, 1 in IDL; 1, 0 in Yorick) I0=1 ; I1=0; // Null value ("" in IDL; [] in Yorick) NULL= []; ioffset= 0; iskip= 0; ncluster= 1; if (param_set(byte_offset)) ioffset= max([0,byte_offset]); if (param_set(byte_skip)) iskip= max([0,byte_skip]); if (param_set(rec_cluster)) ncluster= max([1,rec_cluster]); input_flag= param_set(input); if (typeof(fname_or_handle) == "string") { // Open file if (input_flag) { fhandle= open(fname_or_handle, "rb") //IDL2YORICK: openr, fhandle, fname_or_handle, get_lun=1 ; } else { fhandle= open(fname_or_handle, "w+b") //IDL2YORICK: openw, fhandle, fname_or_handle, get_lun=1 ; } // Return file handle fname_or_handle= fhandle; } else { fhandle= fname_or_handle; } // File size in bytes fsize= sizeof(fhandle); // Data bytes per word bytes_per_word= sizeof(arr(I0)); // Data dimensions ddims= dimsof(arr); ndim= ddims(I0); // Total no. of data words and number of records if (ndim == 0) { nword= 1; nrec= 1; } else { nword= prod(ddims(I0+1:)); if (iskip > 0) nrec= ddims(I0+ndim) ; else nrec= 1; } // Record length (in words/bytes) rec_words= nword/nrec; rec_bytes= rec_words*bytes_per_word; //IDLbegin: //:if ((iskip % rec_bytes) != 0) //: error, "Only complete records may be skipped in IDL"; // No. of records to be skipped //:nskip= iskip/rec_bytes; // Associate record array with file //:frec= assoc( fhandle, replicate(arr(I0), rec_words), ioffset ); // Read/write records //:for (irec=0; irec <= nrec-1; irec++) { //: jrec= irec + (irec/ncluster)*nskip; //: k0= irec*rec_words; //: k1= k0 + rec_words-1; //: if (input_flag) { //: if ((jrec+1)*rec_bytes > (fsize-ioffset)) { // End-of-file reached //: arr= NULL; //: } else { // Read record //: arr(k0:k1)= frec(jrec); //: } //: } else { // Write record //: frec(jrec)= arr(k0:k1); //: } //:} //IDLend: //YORICKbegin: if (nrec == 1) { // Read/write single record if (input_flag) { if ((ioffset+sizeof(arr)) > fsize) arr= NULL; else _read, fhandle, ioffset, arr; } else { _write, fhandle, ioffset, arr; } } else { // Read/write multiple records recarr= array(structof(arr), rec_words); for (irec=0; irec<=nrec-1; irec++) { iaddress= ioffset + irec*rec_bytes + (irec/ncluster)*iskip; if (input_flag) { // Read record if ((iaddress+sizeof(recarr)) > fsize) arr= NULL; else { _read, fhandle, iaddress, recarr; arr(*,I0+irec)= recarr; } } else { // Write record _write, fhandle, iaddress, arr(*,I0+irec); } } } //YORICKend: if (param_set(close_file)) { // Close file close, fhandle; //IDL2YORICK: free_lun, fhandle ; } return; } func oscommand( command, filename=, noprint=) /* DOCUMENT oscommand, command, filename=, noprint=0/1 * Execute the operating system command(s) contained in string COMMAND. * If FILENAME is specified, substitute any % characters in COMMAND with * filename header (i.e., FILENAME minus any suffix and/or leading pathnames). * If FILENAME is an array and the % character occurs in COMMAND, * execute COMMAND once after substituting for each element of the array. * If NOPRINT != 1, suppress printing of the command before execution. * SEE ALSO: */ { func_name= "oscommand"; timer_call,func_name; // Compatibility parameters (including error handling) // compatible.pro: // Parameters for IDL/Yorick compatibility (including error handling) // Array starting index offset, and final index negative offset // (0, 1 in IDL; 1, 0 in Yorick) I0=1 ; I1=0; // Null value ("" in IDL; [] in Yorick) NULL= []; if (is_null(filename) || (strpos(command,"%") < 0)) { // No file name substitution if (!param_set(noprint)) write, command; system, command; } else { nfiles= numberof(filename); for (j=I0; j <= nfiles-I1; j++) { // Extract file name header chars= strsplit( filename(j), "" ); where_slash= where(chars == "/"); if (is_where(where_slash)) chars= chars(where_slash(numberof(where_slash)-I1)+1:); where_dot= where(chars == "."); if (is_where(where_dot)) chars= chars(I0:where_dot(numberof(where_dot)-I1)-1); fheader= strcombine( chars, "" ); // Substitute file name in command chars= strsplit( command, "" ); where_sub= where(chars == "%"); if (is_where(where_sub)) chars(where_sub)= fheader; new_command= strcombine( chars, "" ); if (!param_set(noprint)) write, new_command; system, new_command; } } return timer_return(func_name); } // Maximum calling stack size (for timing info) TIMER_MAXCALL= 256; TIMER_STACK= array(long,TIMER_MAXCALL); TIMER_START= array(double,3,TIMER_MAXCALL); TIMER_NEST= array(double,3,TIMER_MAXCALL); // Maximum recursion level (for timing info) TIMER_RECLEV= 10; TIMER_NAMES= []; struct timer_struc{ // Timer structure long count, curlev; double accum(3,TIMER_RECLEV); } func timer_init(track=,ignore=) /* DOCUMENT timer_init,track=,ignore= Initialize (or reset) execution time counters for functions. Should be called from the main program. Calls to function name strings TRACK, if specified, are stratified based upon the calling function. Calls to function name strings IGNORE, if specified, are ignored (i.e., their timings are lumped with that of the calling function) SEE ALSO: timer_call, timer_return, timer_show, timer */ { extern TIMER_STACK, TIMER_START, TIMER_NEST, TIMER_TRACK, TIMER_IGNORE, TIMER_CALLNO, TIMER_CURFUN, TIMER_CALLFUN, TIMER_NAMES, TIMER_PARS; // Initialize calling stack index TIMER_CALLNO= 0; // Initialize timing operations TIMER_NAMES= "*main*"; TIMER_PARS= timer_struc(); TIMER_TRACK= track; TIMER_IGNORE= ignore; TIMER_CURFUN= ""; TIMER_CALLFUN= ""; timer_call, "*main*"; return; } func timer_call( func_name ) /* DOCUMENT timer_call( func_name ) Log a function call to FUNC_NAME for execution timing. (If FUNC_NAME == "", ignore this call) SEE ALSO: timer_init, timer_return, timer_show, timer */ { extern TIMER_STACK, TIMER_START, TIMER_NEST, TIMER_TRACK, TIMER_IGNORE, TIMER_CALLNO, TIMER_CURFUN, TIMER_CALLFUN, TIMER_NAMES, TIMER_PARS; if (is_void(TIMER_NAMES) || (func_name == "") || anyof(TIMER_IGNORE == func_name) ) return; if (anyof(TIMER_TRACK == func_name)) { stack_name= TIMER_CALLFUN+"->"+func_name; if (stack_name != TIMER_CURFUN) stack_name= TIMER_CURFUN+"->"+func_name; } else { stack_name= func_name; } // Set name of calling function (excluding immediate recursive calls) if (TIMER_CURFUN != stack_name) TIMER_CALLFUN= TIMER_CURFUN; // Name of current function call TIMER_CURFUN= stack_name; // Locate function name list= where(TIMER_NAMES == stack_name); if (!numberof(list)) { // Add new function name to list grow, TIMER_NAMES, stack_name; grow, TIMER_PARS, timer_struc(); ifun= numberof(TIMER_NAMES); } else { ifun= list(1); } // Stack call TIMER_CALLNO += 1; if (TIMER_CALLNO > TIMER_MAXCALL) error, "Too many function calls"; TIMER_STACK(TIMER_CALLNO)= ifun; // Save start time for this call elapsed= array(double,3); timer, elapsed; TIMER_START(,TIMER_CALLNO)= elapsed; TIMER_NEST(,TIMER_CALLNO)= 0.; // Increment calling level and count for function curlev= TIMER_PARS(ifun).curlev; curlev+= 1; if (curlev > TIMER_RECLEV) error, "Too many recursion levels for function '"+func_name+"'"; TIMER_PARS(ifun).curlev= curlev; TIMER_PARS(ifun).count= TIMER_PARS(ifun).count + 1; return; } func timer_return( func_name, return_value ) /* DOCUMENT timer_return( func_name, return_value ) Log a return from function FUNC_NAME for execution timing, with optional RETURN_VALUE, which is simply returned as it is. (If FUNC_NAME == "", ignore this call, and return RETURN_VALUE) (May also be invoked as a subroutine) SEE ALSO: timer_init, timer_call, timer_show, timer */ { extern TIMER_STACK, TIMER_START, TIMER_NEST, TIMER_TRACK, TIMER_IGNORE, TIMER_CALLNO, TIMER_CURFUN, TIMER_CALLFUN, TIMER_NAMES, TIMER_PARS; if (is_void(TIMER_NAMES) || (func_name == "") || anyof(TIMER_IGNORE == func_name) ) return return_value; if (anyof(TIMER_TRACK == func_name)) stack_name= TIMER_CALLFUN+"->"+func_name; else stack_name= func_name; ifun= TIMER_STACK(TIMER_CALLNO); if (TIMER_NAMES(ifun) != stack_name) error, "Called '"+TIMER_NAMES(ifun)+"' but returning from '"+stack_name; // Compute time spent in this function elapsed= TIMER_START(,TIMER_CALLNO); split= array(double,3); timer, elapsed, split; // Accumulate time spent in this function (but not in nested calls) curlev= TIMER_PARS(ifun).curlev; if (curlev < 1) error, "Attempt to return without calling from function '"+func_name+"'"; TIMER_PARS(ifun).accum(,curlev)= TIMER_PARS(ifun).accum(,curlev) - TIMER_NEST(,TIMER_CALLNO) + split; curlev -= 1; TIMER_PARS(ifun).curlev= curlev; // Pop call, and accumulate time spent on nested call for calling function TIMER_CALLNO -= 1; TIMER_NEST(,TIMER_CALLNO)= TIMER_NEST(,TIMER_CALLNO) + split; // Reset name of current function call TIMER_CURFUN= TIMER_NAMES(TIMER_STACK(TIMER_CALLNO)); // Set name of calling function (excluding immediate recursive calls) icallno= TIMER_CALLNO-1; while ((icallno > 1) && (TIMER_CURFUN==TIMER_NAMES(TIMER_STACK(icallno)))) icallno--; TIMER_CALLFUN= TIMER_NAMES(TIMER_STACK(icallno)); return return_value; } func timer_show( ndisp, recursion=, sort_count= ) /* DOCUMENT timer_show, ndisp, recursion=0/1, sortcount=0/1 Display sorted timing statistics for all logged function calls, starting from the function that takes the most execution time. The optional parameter NDISP determines the number of functions to be displayed (the default is to display the top 20). If RECURSION==1, then execution times for different levels of recursive function calls are also reported separately (in addition to the total for the function). If SORT_COUNT==1, then sort by number of calls. NOTE: Execution times for functions whose calls are not logged using TIMER_CALL/TIMER_RETURN are lumped together with the execution times for the functions within which these calls originated. For calls from the main program, unlogged function call timings are lumped together with the execution time for *main*. SEE ALSO: timer_init, timer_call, timer_return, timer */ { extern TIMER_STACK, TIMER_START, TIMER_NEST, TIMER_TRACK, TIMER_IGNORE, TIMER_CALLNO, TIMER_CURFUN, TIMER_CALLFUN, TIMER_NAMES, TIMER_PARS; if (is_void(ndisp)) ndisp= 20; if (is_void(TIMER_NAMES)) error, "Timing counters should be initialized by calling TIMER_INIT"; // Time main program elapsed= TIMER_START(,1); split= array(double,3); timer, elapsed, split; // Compute time spent in main program (but not in nested calls) TIMER_PARS(1).accum(,1)= split - TIMER_NEST(,1); countall= TIMER_PARS().count; accumall= TIMER_PARS().accum; if ((dimsof(accumall))(1) == 2) accumall= accumall(..,-); write, format="%29s %10s %10s %10s %10s\n", "Function", "CPU sec", "System sec", "Wall sec", "Call count" write, format="%29s %10.3f %10.3f %10.3f\n", "*TOTAL*", accumall(1,sum,sum), accumall(2,sum,sum), accumall(3,sum,sum) // Sort in reverse order of execution times if (param_set(sort_count)) { isort= sort( -countall ); } else { isort= sort( -accumall(1,sum,) ); } nfun= min([numberof(TIMER_NAMES), ndisp]); for (j=1; j<=nfun; j++) { ifun= isort(j); // Display total time for function calls write, format="%29s %10.3f %10.3f %10.3f %10d\n", TIMER_NAMES(ifun), accumall(1,sum,ifun), accumall(2,sum,ifun), accumall(3,sum,ifun), countall(ifun); nlev= sum( accumall(sum,,ifun) > 0 ); if ((nlev > 1) && param_set(recursion)) { // Display time for different recursion levels for (ilev=1; ilev<=nlev; ilev++) { write, format="%29d %10.3f %10.3f %10.3f\n", ilev, accumall(1,ilev,ifun), accumall(2,ilev,ifun), accumall(3,ilev,ifun); } } } } /* *---------------------------------------------------------------------- * RCS identification: * $Author:$; $Locker:$ * $Revision:$; $Date:$ GMT; $State:$ * $Source:$ *-------------------------------------------------------------------- * For GNU Emacs: * $Msg-digest-checksum: 6fe216e063d43166d5a0881b801a2d70 $ *-------------------------------------------------------------------- * Local Variables: * mode: C * msg-digest-active: t * End: *-------------------------------------------------------------------- */