; docformat = 'rst' ; ; NAME: ; cgDrawVectors ; ; PURPOSE: ; Plots the velocity vectors of particles at their position. Vectors may be ; overplotted onto other IDL graphics plots. ; ;******************************************************************************************; ; ; ; Copyright (c) 2010, by Fanning Software Consulting, Inc. All rights reserved. ; ; ; ; Redistribution and use in source and binary forms, with or without ; ; modification, are permitted provided that the following conditions are met: ; ; ; ; * Redistributions of source code must retain the above copyright ; ; notice, this list of conditions and the following disclaimer. ; ; * Redistributions in binary form must reproduce the above copyright ; ; notice, this list of conditions and the following disclaimer in the ; ; documentation and/or other materials provided with the distribution. ; ; * Neither the name of Fanning Software Consulting, Inc. nor the names of its ; ; contributors may be used to endorse or promote products derived from this ; ; software without specific prior written permission. ; ; ; ; THIS SOFTWARE IS PROVIDED BY FANNING SOFTWARE CONSULTING, INC. ''AS IS'' AND ANY ; ; EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES ; ; OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT ; ; SHALL FANNING SOFTWARE CONSULTING, INC. BE LIABLE FOR ANY DIRECT, INDIRECT, ; ; INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED ; ; TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; ; ; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ; ; ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT ; ; (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS ; ; SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ; ;******************************************************************************************; ; ;+ ; Plots the velocity vectors of particles at their position. Vectors may be ; overplotted onto other IDL graphics plots. ; ; The program requires the `Coyote Library ` ; to be installed on your machine. ; ; :Categories: ; Graphics ; ; :Params: ; velx: in, required, type=integer/float ; An array containing the X component of the particle velocity vector. ; vely: in, required, type=integer/float ; An array containing the Y component of the particle velocity vector. ; posx_: in, required, type=integer/float ; An array containing the X posiiton of the particle velocity vector. The ; shaft end of the arrow vector is positioned here. ; posy_: in, required, type=integer/float ; An array containing the Y posiiton of the particle velocity vector. The ; shaft end of the arrow vector is positioned here. ; ; :Keywords: ; addcmd: in, optional, type=boolean, default=0 ; An alternative way to set the `Window` keyword. ; clip: in, optional, type=boolean, default=0 ; Set the keyword to clip vector output to the clipping rectangle specified ; in the `CRect` keyword. See the documentation for the IDL graphics ; keyword NOCLIP. ; crect: in, optional, type=float ; A four-element array giving the clipping rectangle [xo, yo, x1, y1]. The ; default clipping rectangle is the plot area. See the documenation for the ; IDL graphics keyword CLIP. ; device: in, optional, type=boolean, default=0 ; Set this keyword to indicate the vector positions are given in ; device coordinates. Data coordinates are assumed. ; fraction: in, optional, type=float, default=1.0 ; A number between 0.0 and 1.0 indicating the fraction of the vectors to ; draw on the plot. Vectors are selected randomly from the total population, ; unless the `Ordered` keyword is set, in which case they are selected ; in an ordered, systematic fashion. For example, Fraction=0.5 will select ; every other input vector. ; hsize: in, optional, type=float ; The size of the the arrow head. By default 1/100th the width ; of the device. (!D.X_SIZE / 100.) ; hthick: in, optional, type=float ; The thickness of the line drawing the arrowhead of the vector. The ; default is 3 for the PostScript device and 1 for all other devices. ; length: in, optional, type=float, default=0.075 ; The length of the `ReferenceVector` in normalized units. All vectors ; are scaled according to this length. ; linestyle: in, optional, type=integer, default=0 ; The line style of the arrow. Line style integers as in PLOT. ; mapCoord: in, optional, type=object ; A map coordinate object (e.g., cgMap) that describes the map projection ; and datum used to specify the vector locations. Note that this could also be a ; map structure as returned from MAP_PROJ_INIT, but in that case the user is ; responsible for setting up the XY map coordinate space independently and ; outside of this program. This coordinate object will be used to transform ; lat/lon locations into the XY locations of the map projection. ; normal: in, optional, type=boolean, default=0 ; Set this keyword to indicate the vector positions are given in ; normalized coordinates. Data coordinates are assumed. ; ordered: in, optional, type=boolean, default=0 ; If this keyword is set, and the `Fraction` keyword is used, the fraction ; of vectors used in the plot are chosen from the entire population in a ; systematic, ordered fashion. ; overplot: in, optional, type=boolean, default=0 ; Set this keyword to overplot the vectors on an established coordinate ; system plot. ; referencevector: in, optional, type=float ; The magnitude of a reference vector that is used to scale all other vectors before display. ; If undefined, 75 percent of the maximum vector magnitude. ; solid: in, optional, type=boolean, default=0 ; Set this keyword to draw solid, filled arrows. ; thick: in, optional, type=float ; The thickness of the line drawing the shaft of the arrow. The ; default is 3 for the PostScript device and 1 for all other devices. ; veccolors: in, optional ; A scalar or vector of colors the same size as `velx`. May be bytes, short integers, ; or strings. Bytes and short integers are treated as indices into the current color ; table. The default is "opposite". ; window: in, optional, type=boolean, default=0 ; Set this keyword to add the command to an cgWindow application. ; xrange: in, optional, type=float ; If a plot is to be drawn, the X range of the plot. By default, the X range is ; calculated to contain all of the velocity vectors. ; yrange: in, optional, type=float ; If a plot is to be drawn, the Y range of the plot. By default, the Y range is ; calculated to contain all of the velocity vectors. ; _extra: in, optional ; Any keywords appropriate for the cgPlot command can be used to create the plot. ; ; :Examples: ; Generate some particle positions and velocities:: ; ; posx = RandomU(seed,200) ; posy = RandomU(seed,200) ; velx = RandomU(seed,200)-0.5 ; vely = RandomU(seed,200)-0.5 ; ; Plot the particle velocities:: ; ; cgDrawVectors, velx, vely, posx, posy ; ; Example using vector colors:: ; ; magnitude = SQRT(velx^2 + vely^2) ; cgLoadCT, 5 ; colors = BytScl(magnitude) ; cgDrawVectors, velx, vely, posx, posy, VecColors=colors ; ; Example using a map projection:: ; ; posx = Indgen(10)*36 ; posy = Replicate(30,10) ; velx = Replicate(10.,10) ; Wind out of the West ; vely = Replicate(0, 10) ; No vertical component ; cgDisplay ; mapNorth = Obj_New('cgMap', 'Polar Stereographic', $ ; Limit=[0, -180, 90, 180], /NoBorder) ; mapNorth -> Draw ; cgMap_Continents, Color='black', Map_Structure=mapNorth ; cgMap_Grid, LatDel=15, LonDel=15, LineStyle=1, Color='charcoal', $ ; /Label, LonLabel=2, Map_Structure=mapNorth ; cgDrawVectors, velx, vely, posx, posy, VecColors='blu6', /Overplot, $ ; ReferenceVector=10, /Solid, Thick=2, MapCoord=mapNorth ; ; :Author: ; FANNING SOFTWARE CONSULTING:: ; David W. Fanning ; 1645 Sheely Drive ; Fort Collins, CO 80526 USA ; Phone: 970-221-0438 ; E-mail: david@idlcoyote.com ; Coyote's Guide to IDL Programming: http://www.idlcoyote.com ; ; :History: ; Change History:: ; Written by David Fanning, based on NASA Astronomy Library program PartVelVec, 22 March 2014 ; Added ORDERED keyword. 24 March 2014. DWF. ; Drawing the vectors All Wrong. Using what I think is the correct algorithm now. 24 March 2014. DWF. ; Modified so I don't change into position variables. 24 March 2014. DWF. ; Fixed a small problem with the algorithm for calculating the direction of the scaled vectors and ; added CLIP and CRECT keywords. 25 March 2014. DWF. ; Added MAPCOORD keyword. 27 March 2014. DWF. ; Modified to set color information once here, instead of allowing cgArrow to do it. 27 March 2014. DWF. ; Some map projections (e.g., polar stereographic) can change the direction of the vector on the ; projected map. To solve this problem, I needed to create a second point in the original coordinate ; system and project both points into the map coordinate system before calculating the angle ; between points. Also then had to figure out how to scale the moved point to the reference ; vector. All appears normal now. 2 Nov 2014. DWF. ; Still a couple of problems in the direction of the vector when the scale is different in the X and ; Y directions. This version of the program allows the vector to be distorted by scale and by a map ; projection. It more closely resembles the NASA program PartVelVec now than it did previously. ; 18 Feb 2015. DWF. ; “To kill an error is as good a service as, and sometimes even better than, the establishing of a ; new truth or fact,” asserted Charles Darwin. In that vein, I have completely gutted the internal ; vector placing algorithm in cgDrawVectors and replaced it with an algorithm that essentially ; duplicates the vectors of the NASA Astronomy Library routine PartVelVec. The essential change was ; to take into account the plot scale in both the X and Y direction when computing the end point ; of a vector. You will see an "aspect" variable in the code. When being placed on map projections, ; using either cgMap_Set or cgMap, the aspect ratio of 180./360. is always used. The aspect ratio is ; used to scale the vector components before placement of the vector. 23 February 2015. DWF. ; Well, duh! Had the WINDOW keyword set when I passed this to a cgWindow. This started a beautiful ; infinite loop! Removed. 11 Nov 2016. DWF. ; Fixed a problem with default vector colors not being byte values. 11 Nov 2016. DWF. ; ; :Copyright: ; Copyright (c) 2014-2016, Fanning Software Consulting, Inc. ;- PRO cgDrawVectors, velx, vely, posx_, posy_, $ ADDCMD=addcmd, $ CLIP=clipit, $ CRECT=crect, $ DEVICE=device, $ FRACTION=fraction, $ HSIZE=hsize, $ HTHICK=hthick, $ LENGTH=length, $ LINESTYLE=linestyle, $ MAPCOORD=mapcoord, $ NORMAL=normal, $ ORDERED=ordered, $ OVERPLOT=overplot, $ REFERENCEVECTOR=referenceVector, $ SOLID=solid, $ THICK=thick, $ VECCOLORS=veccolors_in, $ WINDOW=window, $ XRANGE=xrange, $ YRANGE=yrange, $ _EXTRA=extra Compile_Opt idl2 ; Error handling. Catch, theError IF theError NE 0 THEN BEGIN Catch, /CANCEL void = cgErrorMsg() IF N_Elements(currentState) NE 0 THEN cgSetColorState, currentState RETURN ENDIF ; No parameters. Get some help. IF N_Params() EQ 0 THEN BEGIN Print, 'Syntax of cgDrawVectors: cgDrawVectors, velx, vely, posx, posy' RETURN ENDIF ; Need all four parameters. Get some help. IF N_Params() NE 4 THEN BEGIN Print, 'Syntax of cgDrawVectors: cgDrawVectors, velx, vely, posx, posy' Print, 'All four input parameters are required.' RETURN ENDIF ; All input parameters have to be the same length. IF ((N_Elements(velx) + N_Elements(vely) + N_Elements(posx_) + N_Elements(posy_)) / 4.0) NE N_Elements(velx) THEN BEGIN Message, 'Input parameters are not all the same length.' ENDIF ; Do they want this plot in a resizeable graphics window? IF Keyword_Set(addcmd) THEN window = 1 IF Keyword_Set(window) AND ((!D.Flags AND 256) NE 0) THEN BEGIN ; Special treatment for overplotting or adding a command. IF Keyword_Set(overplot) OR Keyword_Set(addcmd) THEN BEGIN cgWindow, 'cgDrawVectors', velx, vely, posx_, posy_, $ CLIP=clipit, $ CRECT=crect, $ DEVICE=device, $ FRACTION=fraction, $ HSIZE=hsize, $ HTHICK=hthick, $ LENGTH=length, $ LINESTYLE=linestyle, $ MAPCOORD=mapcoord, $ NORMAL=normal, $ ORDERED=ordered, $ OVERPLOT=overplot, $ REFERENCEVECTOR=referenceVector, $ SOLID=solid, $ THICK=thick, $ VECCOLORS=veccolors_in, $ XRANGE=xrange, $ YRANGE=yrange, $ ADDCMD=1, $ _STRICT_EXTRA=extra RETURN ENDIF ; Open a new window or replace the current commands, as required. currentWindow = cgQuery(/CURRENT, COUNT=wincnt) IF wincnt EQ 0 THEN replaceCmd = 0 ELSE replaceCmd=1 cgWindow, 'cgDrawVectors', velx, vely, posx_, posy_, $ CLIP=clipit, $ CRECT=crect, $ DEVICE=device, $ FRACTION=fraction, $ HSIZE=hsize, $ HTHICK=hthick, $ LENGTH=length, $ LINESTYLE=linestyle, $ MAPCOORD=mapcoord, $ NORMAL=normal, $ ORDERED=ordered, $ OVERPLOT=overplot, $ REFERENCEVECTOR=referenceVector, $ SOLID=solid, $ THICK=thick, $ VECCOLORS=veccolors_in, $ XRANGE=xrange, $ YRANGE=yrange, $ REPLACECMD=replaceCmd, $ _STRICT_EXTRA=extra RETURN ENDIF ; Doing this in decomposed color. cgSetColorState, 1, Current=thisColorState ; Default keyword values. SetDefaultValue, fraction, 1.0, RANGE=[0.0, 1.0] SetDefaultValue, length, 0.075, RANGE=[0.0, 1.0] SetDefaultValue, overplot, /Boolean SetDefaultValue, solid, /Boolean IF N_Elements(thick) EQ 0 THEN BEGIN thick = (!D.Name EQ 'PS') ? 3 : 1 ENDIF SetDefaultValue, hthick, thick ; Calculated keyword values. magnitudes = SQRT(velx^2 + vely^2) IF N_Elements(referenceVector) EQ 0 THEN referenceVector = Max(magnitudes) * 0.75D referenceVector = Double(referenceVector) ; Do we need a plot? IF ~overplot THEN cgPlot, [1], /NoData, XRange=[Min(posx_), Max(posx_)], $ YRange=[Min(posy_), Max(posy_)], _STRICT_EXTRA=extra ; Color determination is postponed to here, after a plot is drawn. SetDefaultValue, veccolors_in, Byte(cgColor('opposite')) ; Vector colors should be the same length as other vectors. veclength = N_Elements(velx) IF N_Elements(veccolors_in) EQ 1 THEN veccolors = Replicate(veccolors_in, veclength) ELSE veccolors = veccolors_in ; Calculate default head size. IF N_Elements(hsize) EQ 0 THEN hsize = !D.X_SIZE / 100. ; If we have a map coordinate object, convert the lat/lon position coordinates into XY map grid locations. IF N_Elements(mapCoord) THEN BEGIN xy = mapCoord -> Forward(posx_, posy_, /NoForwardFix) px = Reform(xy[0,*]) py = Reform(xy[1,*]) ENDIF ELSE BEGIN px = posx_ py = posy_ ENDELSE ; We need the positions of vectors in a normalized coordinate system CASE 1 OF Keyword_Set(device): BEGIN nc = Convert_Coord(px, py, /Device, /To_Normal) px = nc[0,*] py = nc[1,*] END Keyword_Set(normal): ELSE: BEGIN ; Data coordinates nc = Convert_Coord(px, py, /Data, /To_Normal) px = nc[0,*] py = nc[1,*] END ENDCASE ; We have to take the plot aspect ratio into account when calculating the vector components. IF (!X.Type EQ 3) || (N_Elements(mapCoord) NE 0) THEN BEGIN ; Map projection aspect = 180./360. ENDIF ELSE BEGIN aspect = (!Y.CRange[1] - !Y.CRange[0]) / (!X.CRange[1] - !X.CRange[0]) ENDELSE ; Modify the longest vector component by the aspect ratio of the plot. IF aspect LT 1.0 THEN BEGIN vx = length*(velx/referenceVector) vy = length*(vely/referenceVector/aspect) ENDIF ELSE BEGIN vx = length*(velx/referenceVector*aspect) vy = length*(vely/referenceVector) ENDELSE x1 = px + vx y1 = py + vy ; Are we doing just a fraction of the vectors. IF fraction LT 1.0 THEN BEGIN ; Make sure you can get some. numGood = Long(fraction * vecLength) IF numGood EQ 0 THEN RETURN ; Compute indices of the vectors to plot. IF Keyword_Set(ordered) THEN BEGIN goodIndices = Long(DIndGen(numGood) / (numGood - 1L) * vecLength) ENDIF ELSE BEGIN goodIndices = Long(RandomU(seed, numGood) * vecLength) ENDELSE px = px[goodIndices] py = py[goodIndices] veccolors = veccolors[goodIndices] x1 = x1[goodIndices] y1 = y1[goodIndices] veclength = N_Elements(px) ENDIF ; Draw the vectors. cgArrow, px, py, x1, y1, $ CLIP=clipit, $ COLOR = veccolors, $ CRECT=crect, $ HSIZE = hsize, $ HTHICK = hthick, $ LINESTYLE=linestyle, $ NORMAL = 1, $ SOLID = solid, $ THICK = thick ; Restore color state cgSetColorState, thisColorState END