local flexplot; /* DOCUMENT flexplot.i --- Flexible plotting package List of routines: * flexp: flexible contouring/shading routine * fpanels: set up multiple plot panels per page * fnext: next plot panel * fproject: apply map projection transformation * fcont: low level contouring/shading routine * fvec: vector plotting routine * fpalette: generate color palette * fget_index: get color index from color name/number * fpad: pads regions with undefined (missing) values * flevels: choose "nice" levels for contouring/axis-divisions * * Version: 1.3 * Modified: 16 Nov 1998, R. Saravanan * * SEE ALSO: yodel, fillc */ require, "color.i" require, "style.i" require, "yodel.i" require, "fillc.i" func flexp(field, xcoord, ycoord, title=, subtitle=, xlim=, ylim=, fill=, overlay=, levs=, type=, width=, c_labels=, line_color=, label_color=, label_size=, low_color=, high_color=, demarc=, mix=, csys=, rev=, stack=, cbar=, cb_labels=, scalef=, units=, vformat=, miss_value=, miss_width=, miss_type=, miss_color=, half_miss=, pad=, rotx=, aspect=, proj=, ppars=, pole_fill=, nested=, triangulate=, charsize=, charthick=, charfont=, xtitle=, xtickformat=, xtickv=, xtickname=, ytitle=, ytickformat=, ytickv=, ytickname=, xpars=, ypars=, apars=, toffset=, position=, overplot=, nodata=) /* DOCUMENT flexp(field, xcoord, ycoord, title=, subtitle=, xlim=, ylim=, fill=, overlay=, levs=, type=, width=, c_labels=, line_color=, label_color=, label_size=, low_color=, high_color=, demarc=, mix=, csys=, rev=, stack=, cbar=, cb_labels=, scalef=, units=, vformat=, miss_value=, miss_width=, miss_type=, miss_color=, half_miss=, pad=, rotx=, aspect=, proj=, ppars=, pole_fill=, nested=, triangulate=, charsize=, charthick=, charfont=, xtitle=, xtickformat=, xtickv=, xtickname=, ytitle=, ytickformat=, ytickv=, ytickname=, xpars=, ypars=, apars=, toffset=, position=, overplot=, nodata=) Flexible contouring/shading routine (essentially a combination of FPALETTE and FCONT) Input: field(*,*) -- field values to be plotted xcoord(*), ycoord(*) -- X, Y coordinate values (KEYWORD PARAMETERS: 0=> false (default), 1=> true) title -- plot title subtitle -- plot subtitle xlim -- X coordinate range [xmin, xmax] ylim -- Y coordinate range [xmin, xmax] fill -- if fill != 0, fill contours (<0 for greyshading) overlay -- if set, overlay color filled plots with contours levs -- contour levels array (non-numeric value => default levels) type -- line type array (re-used cyclically) width -- line width array (re-used cyclically) c_labels -- contour label flag array (re-used cyclically) line_color -- color for contour lines (default "black") label_color -- color for contour labels (default "black") label_size -- size for contour labels (default 1.0) low_color -- color for below-range values (string or negative integer) high_color -- color for above-range values (string or negative integer) demarc -- list of demarcation values for the contour data range mix -- list of mixing colors (name strings or negative integers) csys -- mixing color system ("RGB" or "HSV") rev -- if true, the color scale is reversed stack -- if true, the color scale is stacked on top of previous one cbar -- color bar position = [x0, y0, x1, y1] cb_labels -- color bar labels array (flags, re-used cyclically) scalef -- scale factor for dividing data to be plotted (default 1) (if scalef < 0, it is not displayed, and its absolute value is used) units -- units string for field vformat -- format string for printing out field values miss_value -- missing data value (0 => no missing values, the default) miss_width -- line width for demarcating undefined values miss_type -- line type for demarcating undefined values miss_color -- color for missing values (string or negative integer) half_miss -- if defined, shade only half-grid boxes with missing values at all four corners (NOT YET IMPLEMENTED) pad -- if PAD==1, extrapolate nearest neighbour data values into missing regions; if PAD==2, also interpolate to next nearest neighbours; if PAD==(a float/double value), replace nearest neighbour missing values with PAD. rotx -- if defined, rotate X-coordinate by angle ROTX aspect -- aspect ratio for mesh/projection plots (default 1) proj -- projection ("","NHPOLAR","SHPOLAR","MOLLWEIDE",...) (NOTE: projection names may be abbreviated) ppars -- projection parameters ([start_lon, extreme_lat]/...) pole_fill -- if true, fill in polar value for stereographic projection nested -- if true, shade contours in a nested fashion (a bit slower, but reduces the number of shaded polygons) triangulate -- if true, use routine TRICONTOUR, instead of RECCONTOUR (TRICONTOUR is much slower, but much more robust) charsize -- character size (1.0 for default size) charthick -- character thickness (1.0 for default thickness) charfont -- character font name (default "helvetica") athick -- axis thickness xtitle -- X axis legend string ytitle -- Y axis legend string xpars -- IDL parameters [XTICKS, XMINOR, XTYPE, XSTYLE] defaults to [0, 0, 0, 1] ypars -- IDL parameters [YTICKS, YMINOR, YTYPE, YSTYLE] defaults to [0, 0, 0, 1] apars -- IDL parameters [X/YTHICK, TICKLEN] defaults to [1.0, 0.02] toffset -- [ytitle_offset, xtitle_offset, title_offset, stitle_offset] overplot -- if true, suppress frame/panel advancing Output: (Field is contoured/shaded) SEE ALSO: fpalette, fcont, fproject, fpad, flevels, fpanels */ { // Number of predefined palette colors, number of offset colors, // number of levels, color names, color RGB values, and // the entire palette RGB values extern fpalette_npredef, fpalette_noffset, fpalette_nlev, fpalette_colors, fpalette_rgb, fpalette_palette if (is_void(title)) title= ""; if (is_void(proj)) proj= []; if (is_void(miss_value)) miss_value= 0; if (is_void(miss_rect)) miss_rect= 0; if (is_void(nested)) nested= 1; if (is_void(scalef)) scalef= 1.; if (is_void(units)) units= ""; if (is_void(vformat)) vformat= "%f"; if (is_void(mix) && is_void(demarc)) demarc= [0]; if (is_void(cb_labels)) cb_labels= [1, 0]; if (is_void(overplot)) overplot= 0; if (is_void(charsize)) charsize= 1.0; if (is_void(charfont)) charfont= "helvetica"; xpars2= [0, 0, 0, 1]; ypars2= [0, 0, 0, 1]; apars2= [1.0, 0.02]; if (!is_void(xpars)) xpars2(1:numberof(xpars))= xpars(*); if (!is_void(ypars)) ypars2(1:numberof(ypars))= ypars(*); if (!is_void(apars)) apars2(1:numberof(apars))= apars(*); xticks= xpars2(1); xminor= xpars2(2); xtype= xpars2(3); xstyle= xpars2(4); yticks= ypars2(1); yminor= ypars2(2); ytype= ypars2(3); ystyle= ypars2(4); athick= apars2(1); ticklen= apars2(2); if (is_void(toffset)) toffset= [0., 0., 0., 0.]; // If not overplotting, move to next panel/advance frame if (!overplot) fnext; // No data flag nodata= 0; // Fill flag fill2= 0; if (param_set(fill)) fill2= fill; // Determine field dimensions if (is_void(field)) error, "No field supplied"; fdims= dimsof(field); if (fdims(1) == 0) error, "Cannot plot scalar"; fcount= fdims(2:fdims(1)+1); ndims= sum(fcount > 1); if (ndims == 0) error, "Cannot plot scalar"; if (ndims > 2) error, "Field has too many dimensions"; if (ndims == 2) { // 2-D plot fcount= fcount( where(fcount > 1) ); nx= fcount(1); ny= fcount(2); sfield= field; reshape_array, sfield, [2, nx, ny]; // Generate default "logical" coordinates if (is_void(xcoord)) xcoord= span(0., nx-1., nx); if (is_void(ycoord)) ycoord= span(0., ny-1., ny); if ( ((dimsof(xcoord))(1) == 1) && ((dimsof(ycoord))(1) == 1)) { // 1-D coordinate arrays mesh_flag= 0; if ( (nx != numberof(xcoord)) || (ny != numberof(ycoord)) ) error, "Incorrect X/Y coordinate length"; } else { // 2-D coordinate arrays (mesh) mesh_flag= 1; if ((nx*ny != numberof(xcoord))||(nx*ny != numberof(ycoord))) error, "Incorrect X/Y mesh size"; if ( (is_void(xlim)) || (is_void(ylim)) ) error, "xlim, ylim should be defined for mesh plots"; } } else { // 1-D plot; determine X/Y coordinate dimensions nx= numberof(xcoord); ny= numberof(ycoord); mesh_flag= 0; if (nx*ny != numberof(field)) error, "Incorrect X/Y coordinate length"; if (nx == 1) { sfield= field; reshape_array, sfield, [2, 1, ny]; } if (ny == 1) { sfield= field; reshape_array, sfield, [2, nx, 1]; } } // Copy X/Y coordinates xc= xcoord; yc= ycoord; if (param_set(rotx)) { // Rotate X-coordinate xrotate,360.0,rotx,xc,irot; sfield= rangeop(sfield,"rot",1,count=irot); } if (!mesh_flag) { if ( !is_void(xlim) && (nx > 1) ) { // Choose X subrange xindex= rangeloc( xc, xlim ); if (numberof(xindex) != 2 ) error, "Invalid X-limits" ix0= xindex(1); ix1= xindex(2); } else { ix0= 1; ix1= nx; } if ( !is_void(ylim) && (ny > 1) ) { // Choose X subrange yindex= rangeloc( yc, ylim ); if (numberof(yindex) != 2 ) error, "Invalid Y-limits" iy0= yindex(1); iy1= yindex(2); } else { iy0= 1; iy1= ny; } // Restrict X-Y domain of plot nx= ix1-ix0+1; ny= iy1-iy0+1; xc= float( xc(ix0:ix1) ); yc= float( yc(iy0:iy1) ); } // Copy field in restricted domain sfield= sfield(ix0:ix1,iy0:iy1); // Periodicity flag in X direction xperiodic= 0; if (!is_void(proj)) { // Projection specified projlist= ["mollweide", "nhpolar", "shpolar"]; iproj= strloc( projlist, proj, case_fold=1, abbrev=1, comment="proj" ); projname= projlist(iproj); if (projname == "mollweide") { // Mollweide projection; set default longitude origin if (is_void(ppars) && (!is_void(xlim))) ppars= [ 0.5*(xlim(1)+xlim(2)) ]; } else if ((projname == "nhpolar") || (projname == "shpolar")) { // Polar stereographic plot if (mesh_flag) error, "Projection incompatible with mesh"; // Extend X-coordinate periodically xperiodic= 1; xc= grow(xc(), xc(1)+360.); // Copy array temfield= sfield; if (param_set(pole_fill)) { // Fill pole with mean value sfield= array(float,nx+1,ny+1) if (projname == "nhpolar") { // NH polar stereographic sfield(1:nx,1:ny)= temfield(,); // "North pole value" if ( (miss_value == 0.) || ( sfield(1:nx,ny)(max) < miss_value ) ) sfield(1:nx,ny+1)= sum(sfield(1:nx,ny))/nx; else sfield(1:nx,ny+1)= miss_value; // Extend Y-coordinate to include pole yc= grow(yc(), 90.); } else { // SH polar stereographic sfield(1:nx,2:ny+1)= temfield(,); // "South pole value" if ( (miss_value == 0.) || ( sfield(1:nx,2)(max) < miss_value ) ) sfield(1:nx,1)= sum(sfield(1:nx,2))/nx; else sfield(1:nx,1)= miss_value; // Extend Y-coordinate to include pole yc= grow(-90., yc()) } // Impose periodicity in X-direction; update coordinate counts sfield(nx+1,)= sfield(1,) nx= nx+1; ny= ny+1; } else { // Do not fill the pole; impose periodicity in X-direction sfield= array(float,nx+1,ny); sfield(1:nx,1:ny)= temfield(,); sfield(nx+1,)= sfield(1,); nx= nx+1; } temfield= []; } // Generate mesh of mapped coordinates xc= xc(,-:1:ny); yc= yc(-:1:nx,); mesh= fproject( transpose(grow(xc(*,-), yc(*,-))), proj=proj, ppars=ppars); xc(*)= mesh(1,); yc(*)= mesh(2,); mesh= []; mesh_flag= 1; } // Determine coordinate ranges if (!is_void(xlim)) xlim2= xlim; if (!is_void(ylim)) ylim2= ylim; if ((!is_void(proj)) || mesh_flag) { // Mesh/projection xlim2= [min(xc), max(xc)]; ylim2= [min(yc), max(yc)]; // Determine aspect ratio, if not specified if (is_void(aspect)) aspect= abs( (ylim2(2)-ylim2(1))/(xlim2(2) - xlim2(1)) ); } else { // No mesh/projection if (is_void(xlim)) xlim2= [ "e", "e"]; if (is_void(ylim)) ylim2= [ "e", "e"]; } // Determine max/min field values fmin= min(sfield(*)); fmax= max(sfield(*)); if (miss_value != 0.) { all_undefined= fmin >= miss_value; all_defined= fmax < miss_value; } else { all_undefined= 0; all_defined= 1; } // Padding flags for shading/contouring pad_shade= 0; pad_contour= 0; if ((!is_void(pad)) && (!nodata) && (!all_defined)) { pad_contour= (typeof(pad) == "float") || (typeof(pad) == "double"); pad_shade= pad_contour || param_set(pad); } // Absolute value of scale factor absscalef= abs(scalef); if (all_undefined) { // All data undefined nodata = 1; } else { // Not all undefined values; scale minimum value fmin= fmin / absscalef; if (all_defined) { // No undefined values; scale data fmax= fmax / absscalef; sfield= sfield / absscalef; } else { // Some undefined values occur in data; scale data and find maximum value not_missing = where( sfield < miss_value ); sfield(not_missing)= sfield(not_missing)/absscalef; fmax = max( sfield(not_missing) ); } // Mid-range values fmid = 0.5*(fmax + fmin); } // Determine default contour levels based on the max/min for data, if needed if (!is_number(levs)) levs= flevels( fmin, fmax, color=(fill2 > 0) ); // Determine number of contour levels, demarcation values nlev= numberof(levs); ndemarc= numberof(demarc); ncb_labels= numberof(cb_labels); if (nlev < 1) error, "At least one contour level required"; if ( (nx > 1) && (ny > 1) ) { // 2-D plot // Determine contour line style, if not specified if (is_void(type)) { if ((fill2 <= 0) && (ndemarc == 1)) type= 1 + (levs < demarc(1)); else type= array(1,nlev); } // Construct color bar labelling flag array (extending cyclically) colorlabel= cb_labels; icount= ncb_labels; while (icount < nlev) { colorlabel= [colorlabel, cb_labels]; icount= icount + ncb_labels; } if (is_void(subtitle)) { // No subtitle specified // Create subtitle containing contour interval info subtitle1= swrite( format="("+vformat+" to "+vformat, fmin, fmax); if (numberof(levs) > 1) { cintvl= double(levs(2) - levs(1)); if ((abs(cintvl) >= 0.01) && (abs(cintvl) < 1.)) sintvl= strnum( cintvl, format="%6.4f" ); else if ((abs(cintvl) >= 1.) && (abs(cintvl) < 100.)) sintvl= strnum( cintvl, format="%5.2f" ); else sintvl= strnum( cintvl, format="%9.3e" ); subtitle1= subtitle1 + " by " + sintvl + ")" } else { subtitle1= subtitle1 + ")" } if ((scalef != 1.) && (scalef > 0.)) subtitle1= subtitle1 + " x" + strnum(scalef); if (units != "") subtitle1= subtitle1 + " " + units; } else subtitle1= subtitle; write, format="subtitle=%s\n", subtitle1; // Copy field for possible padding sfieldpad= sfield; if (pad_shade) { // Extrapolate data values into undefined regions miss_rect= 1; // Replace perimeter undefined values with averages over nearest neighbours if (pad_contour) { fpad, sfieldpad, miss_value, miss_value, pad_value=pad; } else { fpad, sfieldpad, miss_value, miss_value, double_pad=(pad > 1); } if (xperiodic) sfieldpad(nx,)= sfieldpad(1,); } else miss_rect= 0; if (fill2 != 0) { // Generate color palette for filling contours fpalette,levs,demarc=demarc,mix=mix,csys=csys,rev=rev,stack=stack, low_color=low_color,high_color=high_color,miss_color=miss_color; } // Shade/contour padded array fcont, sfieldpad, xc, yc, title=title, subtitle=subtitle1, xlim=xlim2, ylim=ylim2, levs=levs, fill=fill, overlay=overlay, type=type, width=width, c_labels=c_labels, miss_value=miss_value, miss_width=miss_width, miss_type=miss_type, miss_rect=miss_rect, aspect=aspect, xperiodic=xperiodic, nested=nested, triangulate=triangulate, line_color=line_color, label_color=label_color, label_size=label_size, ticklen=ticklen, charsize=charsize, charthick=charthick, charfont=charfont, athick=athick, xtitle=xtitle, xtype=xtype, xstyle=xstyle, ytitle=ytitle, ytype=ytype, ystyle=ystyle, xticks=xticks, xminor=xminor, xtickformat=xtickformat, yticks=yticks, yminor=yminor, ytickformat=ytickformat, toffset=toffset } else { // 1-D plot if (!is_void(subtitle)) write, format="subtitle=%s\n", subtitle; if (nx == 1) { // One-dimensional vertical plot if (!superimpose) limits, "e", "e", ylim2(1), ylim2(2) plg, yc, sfield(*), marks=0, width=width, type=type, color=line_color } else { // One-dimensional horizontal plot if (!superimpose) limits, xlim2(1), xlim2(2), "e", "e" plg, sfield(*), xc, marks=0, width=width, type=type, color=line_color } // Display titles port= viewport(); if (!is_null(xtitle) && (xtitle != "")) plt, xtitle, port(zcen:1:2)(1), port(3)-0.050-toffset(2), font=charfont, justify="CT", height=charsize*10.0; if (!is_null(ytitle) && (ytitle != "")) plt, ytitle, port(1)-0.050-toffset(1), port(zcen:3:4)(1), font=charfont, justify="CB", height=charsize*10.0, orient=1; if (!is_null(title) && (title != "")) plt, title, port(zcen:1:2)(1), port(4)+0.02+toffset(3), font=charfont, justify="CB", height=charsize*12.0; if (!is_null(subtitle) && (subtitle != "")) plt, subtitle, port(zcen:1:2)(1), port(3)-0.050- (!is_null(xtitle) && (xtitle != ""))*1.5*charsize*10.*0.0013-toffset(4), font=charfont, justify="CT", height=charsize*10.0; gridxy,0,0,base60=((xtickformat=="dlonlabel")||(xtickformat=="dlatlabel"))+ 2*((ytickformat=="dlonlabel")||(ytickformat=="dlatlabel")); } } func fpanels(nx,ny,file=,pwidth=,pheight=, lmargin=,rmargin=,bmargin=,tmargin=, xgap=,ygap=,landscape=,order=, ptitle=,psubtitle=,toffset=, charcolor=,charsize=,charthick=,charfont=,advance=) /* DOCUMENT fpanels,nx,ny,file=,pwidth=,pheight=, * lmargin=,rmargin=,bmargin=,tmargin=, * xgap=,ygap=,landscape=0/1,order=, * ptitle=,psubtitle=,toffset=, * charcolor=,charsize=,chrathick=,charfont=,advance=0/1 * Sets up multiple (NX*NY) panels on a page, with page width PWIDTH * (default 8.5in) page height PHEIGHT (default 11.0in), left/right/bottom/top * margins LMARGIN/RMARGIN/BMARGIN/TMARGIN (default 1.0in), and X/YGAP * between the panels (default 0.0in). * * If FILE is set to a filename, a new plot window is created, and * PostScript file of that name is used to dump hard copy, by calling HCPON. * * Calling FPANELS with no arguments turns off multiple panels and closes * any associated PostScript file. * * The panels are numbered 1,...,NX*NY as follows: * For ORDER > 0, X ordering varies fastest. * ORDER=1 => Left->right, top->bottom (default) * ORDER=2 => Left->right, bottom->top * ORDER=3 => Right->left, top->bottom * ORDER=4 => Right->left, bottom->top * For ORDER < 0, the X/Y ordering is interchanged, with Y varying fastest. * E.g., * ORDER=-1 => Top->bottom, left->right, and so on. * * (Function FNEXT may be used to select a panel and/or advance to a new page.) * * LANDSCAPE option is for landscape orientation. * * PTITLE/PSUBTITLE specify title/subtitle for the whole page, using characters * of color CHARCOLOR, size CHARSIZE (default 1.0), thickness CHARTHICK * (default 1.0), and font CHARFONT (default "helvetica"). * * TOFFSET=[title_offset, subtitle_offset] controls title positioning. * * Specifying ADVANCE==1 creates a new page for the next plot. * SEE ALSO: fnext, plsys, flexp */ { extern FPANELS_NX, FPANELS_NY, FPANELS_CUR; local landscape1, systems1, legends1, clegends1; if (!param_set(nx) || !param_set(ny)) { // Turn off multiple panels FPANELS_NX= 0; FPANELS_NY= 0; FPANELS_CUR= 0; fma; hcpoff; return; } FPANELS_NX= nx; FPANELS_NY= ny; FPANELS_CUR= 0; // Page width/height (inches) if (is_void(pwidth)) pwidth= 8.5; if (is_void(pheight)) pheight= 11.0; // Left, right, bottom, top margins (inches) if (is_void(lmargin)) lmargin= 1.0; if (is_void(rmargin)) rmargin= 1.0; if (is_void(bmargin)) bmargin= 1.0; if (is_void(tmargin)) tmargin= 1.0; // X/Y gaps if (is_void(xgap)) xgap= 0.0; if (is_void(ygap)) ygap= 0.0; if (is_void(charsize)) charsize= 1.0; if (is_void(charthick)) charthick= 1.0; if (is_void(charfont)) charfont= "helvetica"; if (is_void(toffset)) toffset= [0., 0.]; order1= 1; if (param_set(order)) order1= order; in2ndc = 72.27 * .0013; // convert from inches to NDC coordinate // If hard copy is requested, explicitly create window of right dimensions if (!is_void(file)) { if (strlen(get_env("DISPLAY")) == 0) { display= ""; } else { fma; winkill, 0; display= []; } if (param_set(landscape)) { window,0,width=550,height=425,legends=0,hcp=file,dump=1,display=display; } else { window,0,width=425,height=550,legends=0,hcp=file,dump=1,display=display; } hcpon; } // Read style file read_style, "boxed.gs", landscape1, systems1, legends1, clegends1; // Landscape mode landscape1= int(param_set(landscape)); // Set up default style before copying and changing viewports legends1.nlines = 0; // turn off legends ; clegends1.nlines = 0; // turn off legends ; // Inner dimensions of plot if (landscape1) { // Landscape orientation xinner= pheight - (bmargin+tmargin); yinner= pwidth - (lmargin+rmargin); xoffset= bmargin; yoffset= lmargin; } else { // Portrait orientation xinner= pwidth - (lmargin+rmargin); yinner= pheight - (bmargin+tmargin); xoffset= lmargin; yoffset= bmargin; } // Panel spread xspread= (xinner+0.5*xgap) / nx; yspread= (yinner+0.5*ygap) / ny; // Define viewports systems1= systems1(-:1:nx*ny); xfirst= (order1 > 0); xrightward= (abs(order1) <= 2); yupward= ((abs(order1) % 2) == 0); for (iy=1; iy <= ny; iy++) { for (ix=1; ix <= nx; ix++) { x0= xoffset + (ix-1)*xspread; x1= x0 + xspread - 0.5*xgap; y0= yoffset + (iy-1)*yspread; y1= y0 + yspread - 0.5*ygap; ix1= (2*xrightward-1)*(ix-1 - (1-xrightward)*(nx-1)); iy1= (2*yupward-1)*(iy-1 - (1-yupward)*(ny-1)); ipanel= 1 + (1+xfirst*(nx-1))*iy1 + (1+(1-xfirst)*(ny-1))*ix1; systems1(ipanel).viewport= in2ndc*[x0, x1, y0, y1]; } } // Set style set_style, landscape1, systems1, legends1, clegends1; // Create new page, if requested; if (param_set(advance)) fma; // Display page titles plsys, 0; port= viewport(); if (!is_null(ptitle) && (ptitle != "")) plt, ptitle, port(zcen:1:2)(1), port(4)+0.02+toffset(1), color=charcolor, justify="CB", height=charsize*14.0, font=charfont; if (!is_null(psubtitle) && (psubtitle != "")) plt, psubtitle, port(zcen:1:2)(1), port(3)-0.050-toffset(2), color=charcolor, justify="CT", height=charsize*12.0, font=charfont; return; } func fnext(ipanel,advance=) /* DOCUMENT fnext, ipanel, advance=0/1 * Move to plot panel IPANEL among panels defined by a call to FPANELS. * If none were defined, then simply move to next plot frame. * (If IPANEL > number_of_panels, then 1+(IPANEL-1) % number_panels is used.) * If IPANEL is omitted, the next panel is picked cyclically, and if the * previous panel was the last panel on the page, then a new page is created. * Specifying ADVANCE=1 advances to next page before moving to panel IPANEL, * and if IPANEL is not specified, moves to panel 1 on next page. * SEE ALSO: fpanels, plsys, flexp */ { extern FPANELS_NX, FPANELS_NY, FPANELS_CUR; if (!param_set(FPANELS_NX)) { // Panels not defined; advance frame and return fma; return; } if (!is_void(ipanel)) { FPANELS_CUR= (ipanel-1) % (FPANELS_NX*FPANELS_NY); if (param_set(advance)) fma; plsys, FPANELS_CUR+1; } else { if (param_set(advance)) { advflag= 1; FPANELS_CUR= 1; } else { advflag= (FPANELS_CUR == FPANELS_NX*FPANELS_NY); FPANELS_CUR= 1 + FPANELS_CUR % (FPANELS_NX*FPANELS_NY); } if (advflag) fma; plsys, FPANELS_CUR; } return; } func fproject(x,proj=,ppars=) /* DOCUMENT fproject(x,proj=,ppars=) * Applies the map projection transformation PROJ on a * longitude-latitude degree coordinate array X(2,), returning the * transformed (2,) array. * PROJ="NHPOLAR" => Northern hemisphere polar stereographic (x,y=-1,...,1) * PROJ="SHPOLAR" => Southern hemisphere polar stereographic (x,y=-1,...,1) * PROJ="MOLLWEIDE" => Mollweide projection (x=-2,...,2,y=-1,...,1) * (NOTE: projection names may be abbreviated) * PARS= [start_longitude, extreme_latitude] are option projection parameters * SEE ALSO: fcont, flexp */ { // 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= []; xdims= dimsof(x); if ( (xdims(1) < 1) || (xdims(2) != 2) ) error, "X not a 2-by-n array" radcnv= pi()/180.; // Copy coordinates xproj= x; // Determine projection type projlist= ["mollweide", "nhpolar", "shpolar"]; iproj= strloc( projlist, proj, case_fold=1, abbrev=1, comment="proj" ); projname= projlist(iproj-I1); if (projname == "mollweide") { // Mollweide projection // Default parameter values lon0= 0.; if (numberof(ppars) >= 1) lon0= ppars(I0); xproj(I0,)= cos(x(I0+1,)*radcnv)*(x(I0,)-lon0)/90.; xproj(I0+1,)= sin(x(I0+1,)*radcnv); } else if ((projname == "nhpolar") || (projname == "shpolar")) { // Polar stereographic mapping // Default parameter values lon0= 0.; lat0= 0.; if (numberof(pars) >= 1) lon0= pars(1); if (numberof(pars) >= 2) lat0= pars(2); if (projname == "nhpolar") { // NH polar stereographic mapping lat= [x(2,*), lat0](,max); r0= tan(radcnv*0.5*(90.-lat0)); r= tan(radcnv*0.5*(90.-lat)) / r0; theta= radcnv*(x(1,*) - lon0); } else { // SH polar stereographic mapping lat= [x(2,*), lat0](,min); r0= 1./tan(radcnv*0.5*(90.-lat0)); r= (1./tan(radcnv*0.5*(90.-lat))) / r0; theta= -radcnv*(x(1,*) - lon0); } // Mapped coordinates xproj(1,*)= r * cos(theta); xproj(2,*)= r * sin(theta); } return xproj; } func fcont(field, xcoord, ycoord, title=, subtitle=, xlim=, ylim=, levs=, fill=, overlay=, color=, type=, width=, c_labels=, miss_value=, miss_width=, miss_type=, miss_rect=, aspect=, xperiodic=, nested=, triangulate=, line_color=, label_color=, label_size=, ticklen=, maxseg=, maxpoly=, maxnode=, charsize=, charthick=, charfont=, athick=, xtitle=, xtype=, xstyle=, ytitle=, ytype=, ystyle=, xticks=, xminor=, xtickformat=, yticks=, yminor=, ytickformat=, toffset= ) /* DOCUMENT fcont, field, xcoord, ycoord, title=, subtitle=, xlim=, ylim=, levs=, fill=, overlay=, color=, type=, width=, c_labels=, miss_value=, miss_width=, miss_type=, miss_rect=, aspect=, xperiodic=, nested=, triangulate=, line_color=, label_color=, label_size=, ticklen=, maxseg=, maxpoly=, maxnode=, charsize=, charthick=, charfont=, athick=, xtitle=, xtype=, xstyle=, ytitle=, ytype=, ystyle=, xticks=, xminor=, xtickformat=, yticks=, yminor=, ytickformat=, toffset= Triangulating contouring/shading routine Input: field(*,*) -- field values to be plotted xcoord(*), ycoord(*) -- X, Y coordinate values (KEYWORD PARAMETERS) title -- plot title subtitle -- plot subtitle xlim -- X coordinate range [xmin, xmax] ylim -- Y coordinate range [xmin, xmax] levs -- contour levels array fill -- if set, fill contours overlay -- if set, overlay color filled plots with contours color -- color indices corresponding to contours (NLEV+1 colors are usually required, counting both the below-range and the above-range colors, where NLEV is the number of levels, except when miss_value != 0, when NLEV+2 colors are required, the extra (last) color being used for missing values; FPALETTE may be invoked to generate the list of colors, in which case this parameter should be omitted) type -- contour line type array (re-used cyclically) width -- contour line width array (re-used cyclically) c_labels -- contour label flag array (re-used cyclically miss_value -- missing data value (0 => no missing values, the default) miss_width -- line width for demarcating undefined values miss_type -- line type for demarcating undefined values miss_rect -- if true, mark missing data values using rectangles, rather than using triangles aspect -- aspect ratio for mesh plots (default 1) xperiodic -- if true, data is assumed to be periodic in X direction nested -- if true, shade contours in a nested fashion (a bit slower, but reduces the number of shaded polygons) triangulate -- if true, use routine TRICONTOUR, instead of RECCONTOUR (TRICONTOUR is much slower, but much more robust) line_color -- color for contour lines (default "black") label_color -- color for contour labels (default "black") label_size -- size for contour labels (default 1.0) ticklen -- X/Y axis tick length maxseg -- maximum number of contour segments (default 16384) maxpoly -- maximum number of contour polygons (default 16384) maxnode -- maximum number of contour nodes (default 262144) charsize -- character size (1.0 for default size) charthick -- character thickness (1.0 for default thickness) charfont -- character font name (default "helvetica") athick -- axis thickness xtitle -- X axis legend string xtype -- logarithmic X axis flag ytitle -- Y axis legend string ytype -- logarithmic Y axis flag toffset -- [ytitle_offset, xtitle_offset, title_offset, stitle_offset] Output: (Field is contoured/shaded) SEE ALSO: fproject, fpalette, flevels */ { // Number of levels, number of predefined palette colors, // their names and RGB values, and the entire palette RGB values // Number of predefined palette colors, number of offset colors, // number of levels, color names, color RGB values, and // the entire palette RGB values extern fpalette_npredef, fpalette_noffset, fpalette_nlev, fpalette_colors, fpalette_rgb, fpalette_palette if (is_void(title)) title= ""; if (is_void(type)) type= 1; if (is_void(width)) width= 1.; if (is_void(c_labels)) c_labels= [0, 1]; if (is_void(miss_value)) miss_value= 0; if (is_void(miss_width)) miss_width= 0.; if (is_void(miss_type)) miss_type= 1; if (is_void(miss_rect)) miss_rect= 0; if (is_void(aspect)) aspect= 1.; if (is_void(xperiodic)) xperiodic= 0; if (is_void(nested)) nested= 1; if (is_void(triangulate)) triangulate= 0; if (is_void(line_color)) line_color= -2; if (is_void(label_color)) label_color= -2; if (is_void(label_size)) label_size= 1.0; if (is_void(ticklen)) ticklen= 0.02; if (is_void(charsize)) charsize= 1.0; if (is_void(charfont)) charfont= "helvetica"; if (is_void(maxseg)) maxseg= 16384; if (is_void(maxpoly)) maxpoly= 16384; if (is_void(maxnode)) maxnode= 262144; if (is_void(toffset)) toffset= [0., 0., 0., 0.]; // Logarithmic axes flag logxy, param_set(xtype), param_set(ytype) // No. of line types, line widths, labels ntype= numberof(type); nwidth= numberof(width); nlabels= numberof(c_labels); // Mesh flag mesh_flag= ((dimsof(xcoord))(1) == 2) && ((dimsof(ycoord))(1) == 2); // Plot axes etc. if (mesh_flag) { if ( (is_void(xlim)) || (is_void(ylim)) ) error, "xlim, ylim should be defined for mesh plots"; // Compute aspect ratio of plot area port= viewport(); aspdev= (port(4)-port(3)) / (port(2)-port(1)); // X/Y coordinate multiplication factors if (aspdev <= aspect) { xfac= aspdev / aspect; yfac= 1.; } else { xfac= 1.; yfac= aspect / aspdev; } // Redefine the coordinate system to achieve desired aspect ratio limits, square=0, xlim(1)/xfac, xlim(2)/xfac, ylim(1)/yfac, ylim(2)/yfac } else { // Set plot limits for rectangular plot if (is_void(xlim)) xlim= ["e", "e"]; if (is_void(ylim)) ylim= ["e", "e"]; limits, square=0, xlim(1), xlim(2), ylim(1), ylim(2); } // Determine X/Y coordinate dimensions fdims= dimsof(field); if (fdims(1) != 2) error, "Field not a 2-D array"; nx= fdims(2); ny= fdims(3); // Generate default "logical" coordinates if (is_void(xcoord)) xcoord= span(0., nx-1., nx); if (is_void(ycoord)) ycoord= span(0., ny-1., ny); if (mesh_flag) { if ((nx*ny != numberof(xcoord)) || (nx*ny != numberof(ycoord))) error, "Incorrect X/Y mesh size"; xmesh= float(xcoord); ymesh= float(ycoord); } else { if ( (nx != numberof(xcoord)) || (ny != numberof(ycoord)) ) error, "Incorrect X/Y coordinate length"; // Create mesh coordinate arrays xmesh= float( xcoord(,-:1:ny) ); ymesh= float( ycoord(-:1:nx,) ); } // Determine number of contour levels etc. nlev= numberof(levs); nwidth= numberof(width); nline= numberof(type); if (xperiodic) if (!array_eq(field(1,), field(nx,))) error, "Data not periodic in X direction" // Color values if (param_set(fill)) { if (is_void(color)) { // No colors specified; use colors generated by FPALETTE if ( is_void(fpalette_nlev) || (fpalette_nlev != nlev) ) error, "Incorrect number of contour levels in palette"; // Generate color indices color_ind= char( fpalette_npredef + fpalette_noffset + span(0, nlev+1, nlev+2) ); } else { // Check color specification if (miss_value != 0.) { if (numberof(color) < nlev+2) error, "Not enough colors to deal with missing values"; } else if (numberof(color) < nlev+1) error, "Not enough colors specified"; // Convert color indices to positive "char" values if ( is_void(fpalette_nlev) ) color_ind= color; else color_ind= char(fget_index(color)); } } // Display titles port= viewport(); if (!is_null(xtitle) && (xtitle != "")) plt, xtitle, port(zcen:1:2)(1), port(3)-0.050-toffset(2), font=charfont, justify="CT", height=charsize*10.0; if (!is_null(ytitle) && (ytitle != "")) plt, ytitle, port(1)-0.050-toffset(1), port(zcen:3:4)(1), font=charfont, justify="CB", height=charsize*10.0, orient=1; if (!is_null(title) && (title != "")) plt, title, port(zcen:1:2)(1), port(4)+0.02+toffset(3), font=charfont, justify="CB", height=charsize*12.0; if (!is_null(subtitle) && (subtitle != "")) plt, subtitle, port(zcen:1:2)(1), port(3)-0.050- (!is_null(xtitle) && (xtitle != ""))*1.5*charsize*10.*0.0013-toffset(4), font=charfont, justify="CT", height=charsize*10.0; gridxy,0,0,base60=((xtickformat=="dlonlabel")||(xtickformat=="dlatlabel"))+ 2*((ytickformat=="dlonlabel")||(ytickformat=="dlatlabel")); // Convert variables to proper types for calling C routine if (!triangulate) { // Initialize variables for RECCONTOUR // print, miss_value, maxseg, maxpoly, maxnode polar= 0; // Compute contour fill polygons/contour segments fillp= fillc(field, levs, xmesh, ymesh, contours=!param_set(fill) || param_set(overlay), miss_value=miss_value, miss_rect=miss_rect, xperiodic=xperiodic, nested=nested, maxseg=maxseg, maxpoly=maxpoly, maxnode=maxnode ); if (param_set(fill)) { // Shade polygons; copy vertex coordinates poly_ind= color_ind((*fillp.zp)+1); poly_nodes= *fillp.np; poly_x= *fillp.xp; poly_y= *fillp.yp; // Exclude "transparent" polygons (with color == "bg") cpoly= []; npoly= []; xpoly= []; ypoly= []; icount= 1; for (ipoly=1; ipoly<=numberof(poly_ind); ipoly++) { ncount= poly_nodes(ipoly); if (poly_ind(ipoly) != 0) { grow, cpoly, char(poly_ind(ipoly)); grow, npoly, ncount; grow, xpoly, poly_x(icount:icount+ncount-1); grow, ypoly, poly_y(icount:icount+ncount-1); } icount= icount + ncount; } // Draw filled polygons plfp, cpoly, ypoly, xpoly, npoly; } if ( (!param_set(fill) || param_set(overlay)) && !is_null(fillp.nc) ) { // Draw contours; copy node coordinates vcont= transpose( [*fillp.xc, *fillp.yc] ); // Loop over contours ncnode= 0 for (k=1; k<=numberof(*fillp.zc); k++) { crange= ncnode+1:ncnode+(*fillp.nc)(k); if ((*fillp.zc)(k) < nlev) { // Draw contour plg, vcont(2,crange), vcont(1,crange), marks=0, type=type(1 + (*fillp.zc)(k)%ntype ), width=width(1 + (*fillp.zc)(k)%nwidth ), color=line_color; } else { // Outline missinge value regions // print, k, (*fillp.zc)(k), iseg(k), jseg(k) if (miss_width > 0.) plg, vcont(2,crange), vcont(1,crange), marks=0, type=miss_type, width=miss_width, color=line_color; } // Increment node count ncnode += (*fillp.nc)(k); } } } else { // Initialize variables for TRICONTOUR nvec= 2*(nx-1); maxseg2= (nlev+1)*nvec; maxpoly2= (nlev+1)*nvec; maxnode2= 5*(nlev+1)*nvec; // print, miss_value, maxseg2, maxpoly2, maxnode2; // Initialize variables for TRILIST va= array(float, 3, nvec); vb= array(float, 3, nvec); vc= array(float, 3, nvec); v= array(float, 3, nvec, 3); offset= array(long, 2, 3, 2, 2); // Offset to split along AC diagonal offset(,,,1)= [ [ [0, 0], [1, 1], [0, 1] ], [ [0, 0], [1, 0], [1, 1] ] ]; // Offset to split along BD diagonal offset(,,,2)= [ [ [0, 0], [1, 0], [0, 1] ], [ [1, 0], [1, 1], [0, 1] ] ]; for (j=1; j<=ny-1; j++) { // Shade upper/lower triangles of each row // Initialize vertex coordinates and field values for (i=1; i<=nx-1; i++) { if (miss_value != 0.0) { // Check whether to split rectangle along AC diagonal ac_split= (field(i,j) >= miss_value) || (field(i+1,j+1) >= miss_value); // If miss_rect is non-zero, split rectangle along other diagonal if (miss_rect != 0) ac_split= !ac_split } else { // Split rectangle along AC diagonal ac_split= 1 } for (l=1; l<=3; l++) for (k=1; k<=2; k++) { i2= 2*i-1 + k-1; v(1,i2,l)= xmesh(i + offset(1,l,k,ac_split+1), j + offset(2,l,k,ac_split+1) ); v(2,i2,l)= ymesh(i + offset(1,l,k,ac_split+1), j + offset(2,l,k,ac_split+1) ); v(3,i2,l)= field( i + offset(1,l,k,ac_split+1), j + offset(2,l,k,ac_split+1) ); } } va(,)= v(,,1); vb(,)= v(,,2) vc(,)= v(,,3); // Call TRICONTOUR tricontour, va, vb, vc, nvec, levs, nlev, miss_value, maxseg2, maxpoly2, maxnode2, nseg, lseg, vseg, npoly, ipoly, jpoly, lpoly, nnode, vnode; /* print, "TRICONTOUR: j, nlev, nseg, npoly, nnode:", j, nlev, nseg, npoly, nnode */ vseg= vseg(,,1:nseg); vnode= vnode(,1:nnode); if (param_set(fill)) { // Shade polygons plfp, color_ind(1+lpoly(1:npoly)), vnode(2,1:nnode), vnode(1,1:nnode), jpoly(1:npoly); } if (!param_set(fill) || param_set(overlay)) { // Draw contour segments for (k=1; k<=nseg; k++) { /* print, "VSEG:", k, lseg(k), vseg(1:2,1:2,k) */ if (lseg(k) < nlev) { plg, vseg(2,,k), vseg(1,,k), marks=0, type=type(1+(lseg(k)-1)%ntype), width=width(1+(lseg(k)-1)%nwidth), color=line_color; } else { // Outline undefined values // print, k, lseg(k) if (miss_width > 0.) plg, vseg(2,,k), vseg(1,,k), marks=0, type=miss_type, width=miss_width, color=line_color; } } } } } } func fvec( vy, vx, y, x, scale=, mag=, width=, color=, angle=, head=, arrow=, refvec=) /* DOCUMENT fvec, vy, vx, y, x, scale=, mag=, width=, color=, * angle, head=, arrow=, refvec= * Plot a vector field (VX,VY) on the mesh (X,Y). * The SCALE keyword is the conversion factor from the units of * (VX,VY) to the units of (X,Y) -- a time interval if (VX,VY) is a velocity * and (X,Y) is a position -- which determines the length of the * vector "darts" plotted at the (X,Y) points. If omitted, SCALE is * chosen so that the longest ray arrows have a length comparable * to a "typical" zone size. * The MAG keyword magnifies the arrows (e.g., MAG=2.0 for doubling). * WIDTH controls the line thickness, and COLOR controls the color. * ANGLE is the angle of the arrow head in degrees (defaults to 20 in IDL) * HEAD is the fractional length of the arrow head (defaults to 0.275 in IDL) * ARROW==1 draws Yorick-style arrows (in Yorick only). * REFVEC=[len,x,y] draws a reference arrow of length LEN ending at (X,Y). * SEE ALSO: flexp */ { // 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 (param_set(arrow)) { // Yorick-style arrows //YORICKbegin: if (is_null(angle)) { plv, vy, vx, y, x, scale=scale, width=width, color=color; } else { aspect= tan(angle*(pi()/180.)); plv, vy, vx, y, x, scale=scale, width=width, color=color, aspect=aspect; } if (param_set(mag)) pledit, scalem=mag; return; //YORICKend: } lwidth= 1.; lcolor= 0; if (!is_null(width)) lwidth= width; if (!is_null(color)) lcolor= color; // Angle, fractional length of arrow head anghead= 20. * (pi()/180.); lenhead= .275; if (!is_null(angle)) anghead= angle * (pi()/180.); if (!is_null(head)) lenhead= head; // Number of points npt= numberof(vy); if ( (numberof(vx) != npt) || (numberof(y) != npt) || \ (numberof(x) != npt) ) error, "All arguments should have same dimension"; if (param_set(scale)) { // Specified scale factor scalef= scale; } else { // Determine scale factor ddims= dimsof(vy); nx= ddims(I0+1); ny= ddims(I0+2); // Compute maximum vector length in units of grid interval maxlen= 0.; for (k=I0; k <= ny-I1-1; k++) { for (j=I0; j <= nx-I1-1; j++) { xlen= 0.; ylen= 0.; if (x(j) != x(j+1)) xlen= abs(vx(j,k)/(x(j+1)-x(j))); if (y(k) != y(k+1)) ylen= abs(vy(j,k)/(y(k+1)-y(k))); maxlen= max( [maxlen, sqrt(xlen^2 + ylen^2)] ); } } if (maxlen == 0.) error, "All zero length vectors"; scalef= 1./maxlen; } // Magnify scaling, if requested if (param_set(mag)) scalef= scalef * mag; // Plot vector field on mesh st= lenhead * sin(anghead); ct= lenhead * cos(anghead); for (ipt=I0; ipt <= npt-I1; ipt++) { dx= vx(ipt)*scalef; dy= vy(ipt)*scalef; x0= x(ipt); y0= y(ipt); x1= x0+dx; y1= y0+dy; if ((dx != 0.) || (dy != 0.)) { //IDLbegin: //: plots, [x0, x1, x1-(ct*dx-st*dy), x1, x1-(ct*dx+st*dy)], //: [y0, y1, y1-(ct*dy+st*dx), y1, y1-(ct*dy-st*dx)], //: thick=lwidth, color=lcolor; //IDLend: //YORICKbegin: plg, [y0, y1, y1-(ct*dy+st*dx), y1, y1-(ct*dy-st*dx)], [x0, x1, x1-(ct*dx-st*dy), x1, x1-(ct*dx+st*dy)], marks=0, width=width, color=color //YORICKend: } } if (!is_null(refvec)) { // Plot reference arrow if (numberof(refvec) != 3) error, "Please specify REFVEC=[len,x,y] for reference arrow"; dx= refvec(I0)*scalef; dy= 0.; x1= refvec(I0+1); y1= refvec(I0+2); x0= x1-dx; y0= y1-dy; if ((dx != 0.) || (dy != 0.)) { //IDLbegin: //: plots, [x0, x1, x1-(ct*dx-st*dy), x1, x1-(ct*dx+st*dy)], //: [y0, y1, y1-(ct*dy+st*dx), y1, y1-(ct*dy-st*dx)], //: thick=lwidth, color=lcolor; //IDLend: //YORICKbegin: plg, [y0, y1, y1-(ct*dy+st*dx), y1, y1-(ct*dy-st*dx)], [x0, x1, x1-(ct*dx-st*dy), x1, x1-(ct*dx+st*dy)], marks=0, width=width, color=color //YORICKend: } } return; } func fpalette(levs, demarc=, mix=, csys=, rev=, stack=, low_color=, high_color=, miss_color=) /* DOCUMENT fpalette, levs, demarc=, mix=, csys=, rev=, stack=, low_color=, high_color=, miss_color=) FPALETTE sets up a palette of color indices for filling contour levels. LEVS is the list of contour levels. LOW_COLOR, HIGH_COLOR, and MISS_COLOR are respectively the colors (string or negative integer) to be used for below-range, above-range, and missing data values. If any color is set to "bg", then it is assumed to be transparent. The default values are "bg", "white", and "dark_gray". DEMARC contains a list of demarcation values describing the data range to be contoured. The default is no demarcation values. MIX contains a list of mixing colors (name strings or negative integers). The specified colors are mixed, using the specified color system CSYS ("RGB" or "HSV"). The default value for MIX is ["blue", "white", "red"]. The default value for CSYS is "RGB". if REV==1, then the color scale is reversed. if STACK==1, then the color scale is stacked on top previous one If there are NLEV contour levels, then NLEV-1 mixed colors are generated, one for each contour "band". (NOTE: The last contour level value is actually ignored, because the number of contour "bands" is one less than the number of contour levels.) used to shade the contours. For the special case of a single demarcation value, a list of three mixing colors needs to be specified. E.g., fpalette, [-3,-2,-1,0,1,2], demarc=[0], color=["blue", "white", "red"] associates "white" with the 0 contour, "blue" with negative contours, and "red" with positive contours. If more than one demarcation value is specified, then the number of colors must equal the number of demarcation values, thus establishing a one-to-one correspondence contour values and colors, which is interpolated in a piecewise linear fashion to the generate NLEV-1 colors. There are two additional "predefined" colors available, "light_gray" (-11), and "dark_gray" (-12), in addition to the 10 standard colors. Users may also define their own named colors (numbered -13 onwards) by setting the external variable FPALETTE_EXTRACOL to a list of color names, ant the external variable FPALETTE_EXTRARGB to a 3-by-n array containing RGB triplets (0->255) corresponding to the extra colors. The actual number colors in the generated palette is FPALETTE_NPREDEF+FPALETTE_NOFFSET+NLEV+2, where the external variables FPALETTE_NPREDEF contains the number of predefined colors used by FPALETTE, and FPALETTE_NOFFSET the offset (for stacked color scales) Index FPALETTE_NPREDEF+FPALETTE_NOFFSET+0 corresponds to below range values, index FPALETTE_NPREDEF+FPALETTE_NOFFSET+NLEV to above range values, and index FPALETTE_NPREDEF+FPALETTE_NOFFSET+NLEV+1 to missing values. Indices FPALETTE_NPREDEF+FPALETTE_NOFFSET+1... FPALETTE_NPREDEF+FPALETTE_NOFFSET+NLEV-1 correspond to the first NLEV-1 contour levels. SEE ALSO: flexp, fcont, fget_index */ { // Number of predefined palette colors, number of offset colors, // number of levels, color names, color RGB values, and // the entire palette RGB values extern fpalette_npredef, fpalette_noffset, fpalette_nlev, fpalette_colors, fpalette_rgb, fpalette_palette // Default out of range colors, color mixing scheme if (is_void(low_color)) low_color= "bg"; if (is_void(high_color)) high_color= "white"; if (is_void(miss_color)) miss_color= "dark_gray"; if (is_void(mix)) mix= ["blue", "white", "red"]; if (is_void(csys)) csys= "RGB"; if ( param_set(stack) && !is_void(fpalette_noffset) ) { // No. of offset colors for stacked palettes noffset= fpalette_noffset + fpalette_nlev + 2 } else { noffset= 0; } // Determine number of contour levels, demarcation values, colors nlev= numberof(levs) ndemarc= numberof(demarc); nmix= numberof(mix); if (nlev < 2) error, "At least two contour levels required for palette"; else if (nlev+noffset > 200) error, "Too many contour levels"; if (nmix < 2) error, "At least two mixing colors required for palette"; // Predefined palette color names and RGB values fpalette_colors= [ "bg", "fg", "black", "white", "red", "green", "blue", "cyan", "magenta", "yellow", "light_gray", "dark_gray", "lighter_gray", "darker_gray" ]; fpalette_rgb= [[0,0,0], [255,255,255],[0,0,0], [255,255,255],[255,0,0], [0,255,0],[0,0,255], [0,255,255], [255,0,255],[255,255,0], [192,192,192], [128,128,128], [213,213,213], [96,96,96] ]; // Number of extra predefined colors nextracol= numberof(fpalette_extracol); if (nextracol > 0) { if (!array_eq( dimsof(fpalette_extrargb), [2, 3, nextracol]) ) error, "fpalette_extrargb has wrong structure" // Append extra color names and RGB values fpalette_colors= grow( fpalette_colors, fpalette_extracol); fpalette_rgb= grow( fpalette_rgb, fpalette_extrargb); } // Number of predefined palette colors npredef= numberof(fpalette_colors); // No. of levels in palette, no. of offset colors, and no. of predefined colors fpalette_nlev= nlev; fpalette_noffset= noffset; fpalette_npredef= npredef; if (ndemarc == 0) { // No demarcation values specified; generate evenly spaced values nvals= nmix; cvals= span(levs(1), levs(nlev), nvals); } else if (ndemarc == 1) { // Special case: single demarcation value => "floating" range if (nmix != 3) error,"Three-color mix needed for single demarcation value"; if ((demarc(1) >= levs(1)) && (demarc(1) <= levs(nlev))) { // Demarcation value lies within contour range nvals= 3; clim= max( [demarc(1)-levs(1), levs(nlev)-demarc(1)] ); cvals= demarc(1) + [-clim, 0, clim]; } else { // Demarcation value lies outside contour range; use only two colors if (demarc(1) < levs(1)) cvals= [demarc(1), levs(1), levs(nlev)]; else cvals= [levs(1), levs(nlev), demarc(1)]; nvals= 3; } } else { // Multiple demarcation values specified if (nmix != ndemarc) error, "No. of mixing colors differs from no. of demarcation values"; nvals= ndemarc; cvals= demarc; } // Convert mixing color name/number to RGB values rgb_cmix= transpose( fpalette_rgb(,fget_index(mix)+1) ); // RGB value -> radians conversion factor colrad= 0.5*pi()/255.; if (strmatch(csys,"HSV",1)) { // Convert RGB color table to HSV basis (still called "rgb_cmix", though) rgb_cmix= to_hsv(rgb_cmix); // Enforce periodicity of hue for (i=2; i<=nmix; i++) { if ( (rgb_cmix(i,1)-rgb_cmix(i-1,1)) > 180. ) rgb_cmix(i,1)= rgb_cmix(i,1) - 360.; else if ( (rgb_cmix(i,1)-rgb_cmix(i-1,1)) < -180. ) rgb_cmix(i,1)= rgb_cmix(i,1) + 360.; } } else { // Assume RGB basis; compute cosine transform of RGB values rgb_cmix= cos( rgb_cmix * colrad) ; } // Mix nlev-1 colors rgb_lev= array(float, nlev-1, 3); gray_lev= array(float, nlev-1); k= 1; for (j=1; j<=nlev-1; j++) { // Compute colour mixing fraction while ( (k <= nvals-2) && (levs(j) >= cvals(k+1)) ) k++; frac= float(levs(j) - cvals(k)) / float(cvals(k+1) - cvals(k)); frac= (frac >= 0.) ? frac : 0.; frac= (frac <= 1.) ? frac : 1.; // Modify fraction to avoid pure colors frac= (2./20.) + (17./20.)*frac; // Mix red/green/blue values rgb_lev(j,)= (1.0-frac)*rgb_cmix(k,) + frac*rgb_cmix(k+1,); // Define grayscale gray_lev(j)= float(levs(j) - cvals(1)) / float(cvals(nvals) - cvals(1)); gray_lev(j)= (gray_lev(j) >= 0.) ? gray_lev(j) : 0.; gray_lev(j)= (gray_lev(j) <= 1.) ? gray_lev(j) : 1.; } if (param_set(rev)) { // Reverse color scale rgb_lev= rgb_lev(::-1,); gray_lev= gray_lev(::-1); } if (strmatch(csys,"HSV",1)) { // Revert mixed colors from HSV basis to RGB basis rgb_lev= to_rgb(rgb_lev); } else { // RGB basis; invert cosine transform for mixed colors rgb_lev= acos(rgb_lev) / colrad ; } // Invert cosine transform for grayscale gray_lev= acos(gray_lev) / colrad ; // Palette rgb_palette= array(float,npredef+noffset+nlev+2,3); gray_palette= array(float,npredef+noffset+nlev+2); // Enter predefined palette colors rgb_palette(1:npredef,)= transpose( fpalette_rgb(,) ); gray_palette(1:npredef)= 0. if (noffset > 0) { // Enter previous palette rgb_palette(npredef+1:npredef+noffset,)= fpalette_palette(npredef+1:npredef+noffset,1:3) gray_palette(npredef+1:npredef+noffset)= fpalette_palette(npredef+1:npredef+noffset,4) } // Enter out of range colors rgb_palette(npredef+noffset+1,)= fpalette_rgb(,fget_index(low_color)+1); rgb_palette(npredef+noffset+nlev+1,)=fpalette_rgb(,fget_index(high_color)+1); rgb_palette(npredef+noffset+nlev+2,)=fpalette_rgb(,fget_index(miss_color)+1); gray_palette(npredef+noffset+1)= 255. gray_palette(npredef+noffset+nlev+1)= 255. gray_palette(npredef+noffset+nlev+2)= 0. // Copy mixed level colors rgb_palette(npredef+noffset+2:npredef+noffset+nlev,)= rgb_lev(,); gray_palette(npredef+noffset+2:npredef+noffset+nlev)= gray_lev(); // Convert RGB and gray values to character data type fpalette_palette= grow( char(rgb_palette), char(gray_palette) ); // Load palette palette, fpalette_palette(,1), fpalette_palette(,2), fpalette_palette(,3), fpalette_palette(,4); } func fget_index(color) /* DOCUMENT func fget_index(color) returns the index (0->255) of type char associated with COLOR, which may be a string or a negative integer denoting one of the predefined palette colors. SEE ALSO: fpalette */ { // Number of predefined palette colors, number of offset colors, // number of levels, color names, color RGB values, and // the entire palette RGB values extern fpalette_npredef, fpalette_noffset, fpalette_nlev, fpalette_colors, fpalette_rgb, fpalette_palette ncolor= numberof(color); color_ind= array(char, ncolor); for (i=1; i<=ncolor; i++) { if (typeof(color(i)) == "string") { smatch= strmatch( fpalette_colors, color(i), 1); if (sum(smatch) != 1) error, swrite(format="Invalid/ambiguous color %s\n", color(i)); color_ind(i)= char(smatch(mxx)-1); } else if ( (color(i) > 255) || (color(i) < -numberof(fpalette_colors)) ) error, swrite(format="Out of range color %d\n", color(i)); else { if (color(i) < 0) color_ind(i)= char(-color(i)-1); else color_ind(i)= char(color(i)); } } if (ncolor > 1) return color_ind; else return color_ind(1); } func fpad(field,fmid,missing_value,double_pad=,pad_value=) /* DOCUMENT fpad(field,fmid,missing_value,double_pad=,pad_value=) pads border regions in matrix FIELD between defined and undefined values (i.e., values >= MAX_VALUE) with averaged values (for color shading) and other undefined values with FMID. If DOUBLE_PAD is set, nearest neighbours as well as next nearest neighbours are padded. If PAD_VALUE is defined, use it to pad instead of averaged values. SEE ALSO: flexp */ { next_flag= param_set(double_pad); extrap_flag= is_void(pad_value); // Determine matrix dimensions fdims= dimsof(field); nx= fdims(2); ny= fdims(3); // Compute (padded) defined value mask and matrix n_near= array(long,nx); n_next= array(long,nx); n_near_zero= array(long,nx); n_next_zero= array(long,nx); f_near= array(float,nx); f_next= array(float,nx); f_pad= array(float,nx); idef= array(long,nx+2,3); fdef= array(float,nx+2,3); // Compute mask for "lower" row 1 (j-1) and current row 2 (j) idef(2:nx+1,1)= 0; fdef(2:nx+1,1)= 0.; idef(2:nx+1,2)= field(,1) < missing_value; fdef(2:nx+1,2)= idef(2:nx+1,2) * field(,1); for (j=1; j<=ny; j++) { // Compute masks for "upper" row 3 (j+1) if (j < ny) { idef(2:nx+1,3)= field(,j+1) < missing_value; fdef(2:nx+1,3)= idef(2:nx+1,3) * field(,j+1); } else { idef(2:nx+1,3)= 0; fdef(2:nx+1,3)= 0.; } // Deal with nearest neigbours n_near()=idef(1:nx,2) + idef(3:nx+2,2) + idef(2:nx+1,1) + idef(2:nx+1,3); f_near()=fdef(1:nx,2) + fdef(3:nx+2,2) + fdef(2:nx+1,1) + fdef(2:nx+1,3); // True, if no nearest neighbours n_near_zero()= (n_near() == 0); if (extrap_flag) { // Pad with average of nearest neighbours, or zero f_pad()= f_near() / (n_near()+n_near_zero()); } else { // Pad with specified value f_pad()= pad_value; } if (next_flag) { // Deal with next nearest neigbours n_next()=idef(1:nx,1) + idef(3:nx+2,1) + idef(1:nx,3) + idef(3:nx+2,3); f_next()=fdef(1:nx,1) + fdef(3:nx+2,1) + fdef(1:nx,3) + fdef(3:nx+2,3); // True, if no next nearest neighbours n_next_zero()= (n_next() == 0); if (extrap_flag) { // Further padding with average of next nearest neighbours, or zero f_pad()= f_pad() + n_near_zero()* ( f_next() / (n_next()+n_next_zero()) ); } // Double-pad field field(,j)= fdef(2:nx+1,2) + (1-idef(2:nx+1,2)) * (f_pad() + fmid*n_near_zero()*n_next_zero()); } else { // Single-pad field field(,j)= fdef(2:nx+1,2) + (1-idef(2:nx+1,2)) * (f_pad() + fmid*n_near_zero()); } // Exchange masks idef(2:nx+1,1)= idef(2:nx+1,2) fdef(2:nx+1,1)= fdef(2:nx+1,2) idef(2:nx+1,2)= idef(2:nx+1,3) fdef(2:nx+1,2)= fdef(2:nx+1,3) } } func flevels( xmin, xmax, mindiv=, maxdiv=, color=) /* DOCUMENT flevels( xmin, xmax, mindiv=, maxdiv=, color=) FLEVELS extends any given real interval XMIN... XMAX to an esthetically pleasing interval GMIN...GMAX, and divides it into NDIV divisions of size GDEL each, with NDIV >= MINDIV and NDIV <= MAXDIV. FLEVELS returns the array of NDIV+1 data levels in the range GMIN...GMAX. The default values for MINDIV and MAXDIV are 7 and 15 respectively. If COLOR==1, then the default values are 9 and 19 respectively. FLEVELS would normally be used by a graph plotting program which uses the maxima/minima of the data to choose a scale. SEE ALSO: fpalette, fcont */ { if (is_void(color)) color= 0; if (is_void(mindiv)) mindiv= color ? 9 : 7; if (is_void(maxdiv)) maxdiv= color ? 19 : 15; // Data interval xint = xmax - xmin; // Consider effects of finite machine precision dx= (10.0^(-4)) * max( [abs(xmin), abs(xmax)] ); // Special case : xmin (approximately) equals xmax if (xint <= dx) { if (xmin == 0.0) { gdel= 1.0; gmin= 0.0; gmax= 1.0*mindiv; ndiv= mindiv; return span(gmin, gmax, ndiv+1); } else { // Choose nearest lower power of 10 as data interval gdel= 10.0^floor(log10(abs(xmin))); } } else { // Determine size of a division; start with a power of 10 gdel= 10.0^floor(log10(abs(xint))); // print, "gdel=", gdel, mindiv, maxdiv fmindiv= float(mindiv); fmaxdiv= float(maxdiv); // If too few divisions, decrease division size for (;;) { if ((xint/gdel) >= fmindiv) break; gdel= gdel/2.0; if ((xint/gdel) >= fmindiv) break; gdel= gdel/2.5; if ((xint/gdel) >= fmindiv) break; gdel= gdel/2.0; } // If too many divisions, increase division size for (;;) { if((xint/gdel) <= fmaxdiv) break; gdel = gdel*2.0; if((xint/gdel) <= fmaxdiv) break; gdel = gdel*2.5; if((xint/gdel) <= fmaxdiv) break; gdel = gdel*2.0; } } // Determine the minimum/maximum for the new interval (gmin/gmax) gmin= floor(xmin/gdel)*gdel; gmax= ceil(xmax/gdel)*gdel; ndiv= long( floor(0.5+(gmax-gmin)/gdel) ); // Ensure that there are at least mindiv divisions if (ndiv < mindiv) { ndiv= long( mindiv ); gmax= gmin + ndiv*gdel; } // Compute data levels levels= span(gmin, gmax, ndiv+1); // Locate and reset zero values, if any if ( (levels(1) < 0) && (levels(ndiv+1) > 0) ) for (i=1; i<=ndiv+1; i++) if (abs(levels(i)/gdel) < 10.0^(-4)) levels(i)= 0; return levels; } // MODIFIED HELP_WORKER FUNCTION FROM std.i func help_worker /* xxDOCUMENT help_worker (Not for interactive use -- called by help.) */ { /* help_worker task is pushed by help function -- topic and file arguments are left in help_topic and help_file variables */ topic= help_topic; help_topic= []; file= help_file; help_file= []; if (file) { mark= bookmark(file); line= rdline(file); if (typeof(topic)!="struct_definition") { /* non-struct looks for DOCUMENT comment before any blank lines */ // ***BUG FIX*** n= 30; /* read at most 30 lines looking for DOCUMENT comment */ while (strtok(line)(1) && n--) { if (strmatch(line, "/* DOCUMENT")) break; line= rdline(file); } if (strmatch(line, "/* DOCUMENT")) { do { if (strmatch(line, "**** Y_LAUNCH (computed at runtime) ****")) write, " "+Y_LAUNCH; else if (strmatch(line, "**** Y_SITE (computed at runtime) ****")) write, " "+Y_SITE; else write, line; line= rdline(file); if (!line) break; } while (!strmatch(line, "*/")); write, line; } else { write, ""; } } else { /* struct just prints definition */ gotopen= 0; do { if (!gotopen) gotopen= strmatch(line, "{"); write, line; if (gotopen && strmatch(line, "}")) break; } while (line= rdline(file)); } write, "defined at:"+print(mark)(2); } else { write, ""; info, topic; } } /* *---------------------------------------------------------------------- * RCS identification: * $Author:$; $Locker:$ * $Revision:$; $Date:$ GMT; $State:$ * $Source:$ *-------------------------------------------------------------------- * For GNU Emacs: * $Msg-digest-checksum: 4515f3cdf21186b51b3ce4a6a86bb5a1 $ *-------------------------------------------------------------------- * Local Variables: * mode: C * msg-digest-active: t * End: *-------------------------------------------------------------------- */