PRO STR_FIELD1,U,V,X,Y, Missing=Missing, Length=length, Dots=dots, $ Title=title, position=position, noerase=noerase, color=color,$ xtitle=xtitle, ytitle=ytitle, xmargin=xmargin, ymargin=ymargin, $ easy_view=easy_view, noaxis=noaxis ; ;+ ; NAME: ; STR_FIELD ; ; PURPOSE: ; Produce a two-dimensional velocity field and streamlines plot. ; ; See notes at http://www.amara.com/ftpstuff/streamlines1.txt ; and consider using an integration method better than Euler. ; ; A directed arrow is drawn at each point showing the direction and ; magnitude of the field and streamlines connect the midpoints of the ; arrows. These streamlines are calculated by integrating tangents ; to the vector field ; ; CATEGORY: ; Plotting, two-dimensional. ; ; CALLING SEQUENCE: ; STR_FIELD, U, V [, X, Y] ; ; INPUTS: ; U: The X component of the two-dimensional field. ; U must be a two-dimensional array. ; ; V: The Y component of the two dimensional field. Y must have ; the same dimensions as X. The vector at point (i,j) has a ; magnitude of: ; ; (U(i,j)^2 + V(i,j)^2)^0.5 ; ; and a direction of: ; ; ATAN2(V(i,j),U(i,j)). ; ; OPTIONAL INPUT PARAMETERS: ; X: Optional abcissae values. X must be a vector with a length ; equal to the first dimension of U and V. ; ; Y: Optional ordinate values. Y must be a vector with a length ; equal to the first dimension of U and V. ; ; KEYWORD INPUT PARAMETERS: ; MISSING: Missing data value. Vectors with a LENGTH greater ; than MISSING are ignored. ; ; LENGTH: Length factor. The default of 1.0 makes the longest (U,V) ; vector the length of a cell. ; ; DOTS: Set this keyword to 1 to place a dot at each missing point. ; Set this keyword to 0 or omit it to draw nothing for missing ; points. Has effect only if MISSING is specified. ; ; TITLE: A string containing the plot title. ; ; POSITION: A four-element, floating-point vector of normalized ; coordinates for the rectangular plot window. ; This vector has the form [X0, Y0, X1, Y1], where (X0, Y0) ; is the origin, and (X1, Y1) is the upper-right corner. ; ; NOERASE: Set this keyword to inhibit erase before plot. ; ; COLOR: The color index used for the plot. ; ; EASY_VIEW: Set this keyword if you don't wish to have all streamlines ; shown. For example if you have more than 50 or 100 vectors in your ; vector plot, the streamline plotting will add considerable time to ; your plot. ; ; XMARGIN: A two-element vector that specifies the vertical margin between ; the axis and the screen border in character units. ; ; YMARGIN: A two-element vector that specifies in the horzontal margin ; between the map and screen border in character units. ; ; NOAXIS: Set this keyword if you wish to have no axes drawn (i.e. ; if you will be calling this routine on top of another plot. ; ; OUTPUTS: ; None. ; ; COMMON BLOCKS: ; None. ; ; SIDE EFFECTS: ; Plotting on the selected device is performed. System ; variables concerning plotting are changed. ; ; RESTRICTIONS: ; None. ; ; PROCEDURE: ; Straightforward. The system variables !XTITLE, !YTITLE and ; !MTITLE can be set to title the axes. ; ; MODIFICATION HISTORY: ; VEL_FIELD: DMS, RSI, Oct., 1983. ; ; VEL_FIELD: For Sun, DMS, RSI, April, 1989. ; ; VEL_FIELD: Added TITLE, Oct, 1990. ; ; VEL_FIELD: Added POSITION, NOERASE, COLOR, Feb 91, RES. ; ; Created STR_FIELD as modification of VEL_FIELD, ; Nov/Dec 1993, Amara Graps NASA-Ames ;- ; on_error,2 ;Return to caller if an error occurs s = size(u) t = size(v) if s(0) ne 2 then begin baduv: message, 'U and V parameters must be 2D and same size.' endif if total(abs(s(0:2)-t(0:2))) ne 0 then goto,baduv ; if n_params(0) lt 3 then x = findgen(s(1)) else $ if n_elements(x) ne s(1) then begin badxy: message, 'X and Y arrays have incorrect size.' endif if n_params(1) lt 4 then y = findgen(s(2)) else $ if n_elements(y) ne s(2) then goto,badxy ; if n_elements(missing) le 0 then missing = 1.0e30 if n_elements(length) le 0 then length = 1.0 mag = sqrt(u^2+v^2) ;magnitude. stream = mag ;-------- Subscripts of good elements ---------------------- nbad = 0 ;# of missing points if n_elements(missing) gt 0 then begin good = where(mag lt missing) if keyword_set(dots) then bad = where(mag ge missing, nbad) endif else begin good = lindgen(n_elements(mag)) endelse mag = mag(good) ;Discard missing values maxmag = max(mag) ugood = u(good) vgood = v(good) ;Define grid for vector and streamline plotting xleft = min(x) xright = max(x) ybottom = min(y) ytop = max(y) delta_x = abs(xleft - xright) delta_y = abs(ybottom - ytop) nsteps = 300 ;# of integration steps per streamline delta_max = delta_x > delta_y ;choose greater delta h = delta_max / 30. ;step-size xstream = fltarr(n_elements(ugood), nsteps, 2) ystream = fltarr(n_elements(vgood), nsteps, 2) xmid = fltarr(n_elements(ugood)) ymid = fltarr(n_elements(vgood)) sina = length * (xright-xleft)/s(1)/maxmag*ugood ;sin & cosine components. cosa = length * (ytop-ybottom)/s(2)/maxmag*vgood ;Note: x0, x1, y0, y1 are boundaries for determining streamline limits ; ;print, 'xtitle = ', xtitle ;print, 'ytitle = ', ytitle ;-------------- plot to get axes --------------- ;Set up some of the keywords if n_elements(title) le 0 then title = '' if n_elements(xtitle) le 0 then xtitle = 'U(x,y)' if n_elements(ytitle) le 0 then ytitle = 'V(x,y)' if n_elements(xmargin) le 0 then xmargin = !x.margin if n_elements(ymargin) le 0 then ymargin = !y.margin if n_elements(color) eq 0 then color = !p.color if n_elements(position) ne 0 then begin plot,[xleft,xright],[ytop,ybottom],/nodata,/xst,/yst,title=title, $ noerase=noerase, color=color, xtitle=xtitle, ytitle = ytitle,$ xmargin=xmargin, ymargin = ymargin, position = position endif else $ if keyword_set(noaxis) then begin plot,[xleft,xright],[ytop,ybottom],/nodata,$ noerase=noerase, color=color, xmargin=xmargin, $ ymargin = ymargin, xstyle=14, ystyle=14 endif else $ plot,[xleft,xright],[ytop,ybottom],/nodata,/xst,/yst,$ title=title, noerase=noerase, color=color, xtitle=xtitle, $ ytitle=ytitle, xmargin=xmargin, ymargin = ymargin ; ;--- place vectors with arrows -- r = .3 ;len of arrow head angle = 22.5 * !dtor ;Angle of arrowhead st = r * sin(angle) ;sin 22.5 degs * length of head ct = r * cos(angle) ; for i=0,n_elements(good)-1 do begin ;Each point x0 = x(good(i) mod s(1)) ;get coords of start & end dx = sina(i) x1 = x0 + dx xmid(i) = (x1 + x0)/2.0 y0 = y(good(i) / s(1)) dy = cosa(i) y1 = y0 + dy ymid(i) = (y0 + y1)/2.0 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)], $ color=color, thick = 2.0 endfor if nbad gt 0 then $ ;Dots for missing? oplot, x(bad mod s(1)), y(bad / s(1)), psym=3, color=color ;-----------integrate for streamlines---------- ;Set up for showing a select set of streamlines nskip = 1 if keyword_set(easy_view) then begin nshow_streams = 20. ;# of streamlines to show for easy viewing ntotal = n_elements(good) ;# of potential streamlines (1 per vector) ;determine increment value for "for loop" based on nshow_streams & ntotal if ntotal gt nshow_streams then nskip = fix(ntotal/nshow_streams) end ;if easy_view istep = 0 ;step counter for i = 0, n_elements(good) - 1, nskip do begin ;each point (field vector) for j = 0, 1 do begin ;forward and backward streamlines curr_x = xmid(i) curr_y = ymid(i) xstream(i,istep, j) = curr_x ystream(i,istep, j) = curr_y out_of_bounds = 0 ;false while ((istep le nsteps-2) and (out_of_bounds eq 0)) do begin ;Calculate gx, gy (deriv f'ns of x,y) case 1 of istep eq 0: begin ;1st step only gx = sina(i) gy = cosa(i) end ;case istep eq 0 else: begin ;Calc distance of all other vectors from current point tempx = xmid - curr_x tempy = ymid - curr_y distance = sqrt(tempx^2 + tempy^2) ;Sort to get nearest ones neari = sort(distance) ;array indices nearest = distance(neari) ;sorted by distance from curr nearest_gx = sina(neari) ;function in x direction at i nearest_gy = cosa(neari) ;function in y direction at i ;Choose nearest 5 vectors, and use to calculate a weighted ;average. Our weights for the nearest points are inversely ;proportional to their distance from the current vector. ;The 1st vector nearest(0) is the current point, which we don't ;want. ; ;Equation being solved for the weights is: ;k/d1 + k/d2 + k/d3 +k/d4 = 1; ;d1, d2, d3, d4 = nearby vector distances from current vector ;w1 = k/d1, w2 = k/d2, w3 = k/d3, w4 = k/d4 ; w1, w2, w3, w4 = weights for nearby ;vectors, and k = a proportionality constant which we must calculate. d1 = nearest(1) & d2 = nearest(2) & d3 = nearest(3) & d4 = nearest(4) k = (d1*d2*d3*d4) / (d2*d3*d4 + d1*d3*d4 + d1*d2*d4 + d1*d2*d3) w1 = k/float(d1) & w2 = k/float(d2) & w3 = k/float(d3) & w4 = k/float(d4) ;x-component d1x = nearest_gx(1) & d2x = nearest_gx(2) d3x = nearest_gx(3) & d4x = nearest_gx(4) gx = (w1*d1x + w2*d2x + w3*d3x + w4*d4x)/4.0 ;y-component d1y = nearest_gy(1) & d2y = nearest_gy(2) d3y = nearest_gy(3) & d4y = nearest_gy(4) gy = (w1*d1y + w2*d2y + w3*d3y + w4*d4y)/4.0 end ;case else i endcase istep = istep + 1 ;Apply Eulers algorithm to integrate to next (forward or backward) step in ;streamline case j of 0: begin ;forward xstream(i,istep, j) = curr_x + gx * h ystream(i,istep, j) = curr_y + gy * h end 1: begin ;backward xstream(i,istep, j) = curr_x - gx * h ystream(i,istep, j) = curr_y - gy * h end endcase curr_x = xstream(i,istep, j) curr_y = ystream(i,istep, j) ;See if out of bounds to end integration case 1 of curr_x le xleft: out_of_bounds = 1 ;true curr_x ge xright: out_of_bounds = 1 curr_y le ybottom: out_of_bounds = 1 curr_y ge ytop: out_of_bounds = 1 else: out_of_bounds = 0 ;false endcase end ;while istep le nsteps oplot, xstream(i,0:istep, j), ystream(i,0:istep, j) istep = 0 endfor ;j (forward and backward streamline) endfor ;i (each wind vector) end ;**************************************************************** PRO STR_FIELD2,U,V,X,Y, Missing=Missing, Length=length, Dots=dots, $ Title=title, position=position, noerase=noerase, color=color,$ xtitle=xtitle, ytitle=ytitle, xmargin=xmargin, ymargin=ymargin, $ easyview=easyview, noaxis=noaxis, nolabels=nolabels ; ;+ ; NAME: ; STR_FIELD2 ; ; PURPOSE: ; Produce a two-dimensional velocity field plot with streamlines ; by finding the level sets of stream functions. ; ; See notes at http://www.amara.com/ftpstuff/streamlines2.txt ; ; A directed arrow is drawn at each point showing the direction and ; magnitude of the field and streamlines show the field in a mass- ; conservation way. ; ; CATEGORY: ; Plotting, two-dimensional. ; ; CALLING SEQUENCE: ; STR_FIELD, U, V [, X, Y] ; ; INPUTS: ; U: The X component of the two-dimensional field. ; U must be a two-dimensional array. ; ; V: The Y component of the two dimensional field. Y must have ; the same dimensions as X. The vector at point (i,j) has a ; magnitude of: ; ; (U(i,j)^2 + V(i,j)^2)^0.5 ; ; and a direction of: ; ; ATAN2(V(i,j),U(i,j)). ; ; OPTIONAL INPUT PARAMETERS: ; X: Optional abcissae values. X must be a vector with a length ; equal to the first dimension of U and V. ; ; Y: Optional ordinate values. Y must be a vector with a length ; equal to the first dimension of U and V. ; ; KEYWORD INPUT PARAMETERS: ; MISSING: Missing data value. Vectors with a LENGTH greater ; than MISSING are ignored. ; ; LENGTH: Length factor. The default of 1.0 makes the longest (U,V) ; vector the length of a cell. ; ; DOTS: Set this keyword to 1 to place a dot at each missing point. ; Set this keyword to 0 or omit it to draw nothing for missing ; points. Has effect only if MISSING is specified. ; ; TITLE: A string containing the plot title. ; ; POSITION: A four-element, floating-point vector of normalized ; coordinates for the rectangular plot window. ; This vector has the form [X0, Y0, X1, Y1], where (X0, Y0) ; is the origin, and (X1, Y1) is the upper-right corner. ; ; NOERASE: Set this keyword to inhibit erase before plot. ; ; COLOR: The color index used for the plot. ; ; EASYVIEW: Set this keyword if you don't wish to have all vectors ; shown. For example if you have more than ~200 vectors in your ; vector plot, the plotting will add considerable time. ; ; XMARGIN: A two-element vector that specifies the vertical margin between ; the axis and the screen border in character units. ; ; YMARGIN: A two-element vector that specifies in the horzontal margin ; between the map and screen border in character units. ; ; NOAXIS: Set this keyword if you wish to have no axes drawn (i.e. ; if you will be calling this routine on top of another plot or map. ; ; NOLABELS: Set this keyword if you wish to *not* have the contours ; of the streamlines labeled with their values. ; ; OUTPUTS: ; None. ; ; COMMON BLOCKS: ; None. ; ; SIDE EFFECTS: ; Plotting on the selected device is performed. System ; variables concerning plotting are changed. ; ; RESTRICTIONS: ; None. ; ; PROCEDURE: ; Straightforward for vectors. The streamlines are calculated ; using Bob Chatfield's and John Vastano's mass-conservation method. ; (i.e. producing level sets of the stream function) ; ; MODIFICATION HISTORY: ; VEL_FIELD: DMS, RSI, Oct., 1983. ; ; VEL_FIELD: For Sun, DMS, RSI, April, 1989. ; ; VEL_FIELD: Added TITLE, Oct, 1990. ; ; VEL_FIELD: Added POSITION, NOERASE, COLOR, Feb 91, RES. ; ; Created STR_FIELD as modification of VEL_FIELD, ; Nov/Dec 1993, Amara Graps NASA-Ames ;- ; on_error,2 ;Return to caller if an error occurs s = size(u) t = size(v) if s(0) ne 2 then begin baduv: message, 'U and V parameters must be 2D and same size.' endif if total(abs(s(0:2)-t(0:2))) ne 0 then goto,baduv ; if n_params(0) lt 3 then x = findgen(s(1)) else $ if n_elements(x) ne s(1) then begin badxy: message, 'X and Y arrays have incorrect size.' endif if n_params(1) lt 4 then y = findgen(s(2)) else $ if n_elements(y) ne s(2) then goto,badxy ; if n_elements(missing) le 0 then missing = 1.0e30 if n_elements(length) le 0 then length = 1.0 mag = sqrt(u^2+v^2) ;magnitude. ;-------- Subscripts of good elements ---------------------- nbad = 0 ;# of missing points if n_elements(missing) gt 0 then begin good = where(mag lt missing) if keyword_set(dots) then bad = where(mag ge missing, nbad) endif else begin good = lindgen(n_elements(mag)) endelse mag = mag(good) ;Discard missing values maxmag = max(mag) ugood = u(good) vgood = v(good) ;Define grid for vector and streamline plotting xleft = x(0) xright = x(n_elements(x)-1) ybottom = y(0) ytop = y(n_elements(y)-1) range_x = (xright - xleft) range_y = (ytop - ybottom) ;Set grid for stream function nsteps_x = n_elements(x);# of steps in each dimension nsteps_y = n_elements(y);# of steps in each dimension div_x = nsteps_x-1 div_y = nsteps_y-1 delta_x = range_x / float(div_x) delta_y = range_y / float(div_y) stream = fltarr(nsteps_x, nsteps_y) stream(0,0) = 0.0 ;initialize first point xbeg = fltarr(n_elements(ugood)) ybeg = fltarr(n_elements(vgood)) sina = length * (xright-xleft)/s(1)/maxmag*ugood ;sin & cosine components. cosa = length * (ytop-ybottom)/s(2)/maxmag*vgood ;Note: x0, x1, y0, y1 are streamline boundaries ;-------------- plot to get axes --------------- ;Set up some of the keywords if n_elements(title) le 0 then title = '' if n_elements(xtitle) le 0 then xtitle = 'U(x,y)' if n_elements(ytitle) le 0 then ytitle = 'V(x,y)' if n_elements(xmargin) le 0 then xmargin = !x.margin if n_elements(ymargin) le 0 then ymargin = !y.margin if n_elements(color) eq 0 then color = !p.color if n_elements(position) ne 0 then begin plot,[xleft,xright],[ytop,ybottom],/nodata,/xst,/yst,title=title, $ noerase=noerase, color=color, xtitle=xtitle, ytitle = ytitle,$ xmargin=xmargin, ymargin = ymargin, position = position endif else $ if keyword_set(noaxis) then begin plot,[xleft,xright],[ytop,ybottom],/nodata, $ noerase=noerase, color=color, xstyle=5,ystyle=5,$ xmargin=xmargin,ymargin=ymargin endif else $ plot,[xleft,xright],[ytop,ybottom],/nodata,/xst,/yst,$ title=title, noerase=noerase, color=color, xtitle=xtitle, $ ytitle=ytitle, xmargin=xmargin, ymargin = ymargin ;Set up for showing a select set of vectors nskip = 1 if keyword_set(easyview) then begin nshow_vecs = 200. ;# of vectors to show for easy viewing (EMPIRICAL!) ntotal = n_elements(good) ;# of potential vectors ;Determine increment value for "for loop" based on nshow_vecs & ntotal if ntotal gt nshow_vecs then nskip = fix(ntotal/nshow_vecs) end ;if easyview ; ;----- Calculate vectors ------------- r = .3 ;len of arrow head angle = 22.5 * !dtor ;Angle of arrowhead st = r * sin(angle) ;sin 22.5 degs * length of head ct = r * cos(angle) ; ;First calculate so that we have the full vector xbeg, ybeg for i=0,n_elements(good)-1, 1 do begin ;Each point x0 = x(good(i) mod s(1)) ;get coords of start & end dx = sina(i) x1 = x0 + dx xbeg(i) = x0 y0 = y(good(i) / s(1)) dy = cosa(i) y1 = y0 + dy ybeg(i) = y0 endfor ;-------Place vectors with arrows------------------------ for i=0,n_elements(good)-1, nskip do begin ;Each point x0 = x(good(i) mod s(1)) ;get coords of start & end dx = sina(i) x1 = x0 + dx y0 = y(good(i) / s(1)) dy = cosa(i) y1 = y0 + dy 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)], $ color=color, thick = 1.25,noclip=0 endfor if nbad gt 0 then $ ;Dots for missing? oplot, x(bad mod s(1)), y(bad / s(1)), psym=3, color=color ;-----------Create a streamline grid---------- ;CURRENTLY NOT USED-------------------- ;Set up for showing a select set of streamlines nskip = 1 if keyword_set(easyview) then begin nshow_streams = 20. ;# of vectors to show for easy viewing ntotal = n_elements(good) ;# of potential streamlines ;Determine increment value for "for loop" based on nshow_streams & ntotal if ntotal gt nshow_streams then nskip = fix(ntotal/nshow_streams) end ;if easyview ;-------------------------------------- ;Fill up the Stream function grid by moving "UP" (i.e. constant x) yrow = 1 while yrow le nsteps_y-1 do begin ;Get u{0,n} and u{0,n+1} gx = fltarr(2) gx(0) = u(0,yrow-1) gx(1) = u(0,yrow) ;Use Trapezoidal Rule with the two "gx" evaluations stream(0,yrow) = -0.5*(gx(0) + gx(1))*delta_y + stream(0,yrow-1) yrow = yrow + 1 end ;while filling "UP" ;Now fill up the Stream grid by moving "ROW by ROW" (i.e. constant y) yrow = 0 while yrow le nsteps_y-1 do begin istep = 1 while istep le nsteps_x-1 do begin ;Find v{0,n} and v{0,n+1} gy = fltarr(2) gy(0) = v(istep-1, yrow) gy(1) = v(istep,yrow) ;Use Trapezoidal Rule with the two "gy" evaluations stream(istep,yrow) = 0.5*(gy(0) + gy(1))*delta_x + stream(istep-1,yrow) istep = istep + 1 end; while istep yrow = yrow + 1 end ;while filling "ROW by ROW" (yrow) ;-----------Contour the streamline grid---------- ;Set up levels and colors for contour ncontour = 25 nclevel = ncontour+1 datamin = min(stream) datamax = max(stream) dinc = (datamax-datamin)/ncontour levels = datamin+dinc*findgen(nclevel) clabels = intarr(nclevel) for i=0,ncontour do clabels(i) = 0 for i = 0,ncontour,2 do clabels(i) = 1 cinc = 200./ncontour colors=30.+cinc*findgen(nclevel) ;Now contour if keyword_set(nolabels) then begin contour,stream,x,y,levels=levels, c_colors=colors,/noerase, $ xstyle=5, ystyle=5, xmargin = xmargin, ymargin= ymargin endif else $ contour,stream,x,y,levels=levels, c_colors=colors,c_charsize=1.0,$ c_labels=clabels,/follow, /noerase, xstyle=5, ystyle=5, $ xmargin = xmargin, ymargin= ymargin end ;of procedure str_field2 ;**************************************************************** PRO TEST_STRFIELD ;This program initializes the "u" and "v" ;vector components of a vector field to test the stream ;and velocity field routine; str_field.pro. ;----------------------------------------------------- ;Amara Graps 11-19-93 Updated: 12-30-93 ;---------------------------------------------------- u = fltarr(6,6) v = fltarr(6,6) ;Equation for 2-D vector field is Fvec = -y i + x j ;where i and j are unit vectors pointing in ijk coords ;F1(x,y)i = -y --> our "u" vector ;F2(x,y)j = x --> our "v" vector x=indgen(6) y=indgen(6) ;fill in edge points for i = 0, 5 do u(0,i) = float(-i) for i = 0, 5 do v(i,0) = float(i) ;fill in some intermediate points u(1,1) = -1. for i = 0,5 do u(i,1) = -1. u(2,2) = -2. for i = 0,5 do u(i,2) = -2. u(2,3) = -3. for i = 0,5 do u(i,3) = -3. u(4,4) = -4. for i = 0,5 do u(i,4) = -4. u(5,5) = -5. for i = 0,5 do u(i,5) = -5. v(1,1) = 1. for i = 0,5 do v(1,i) = 1. v(2,2) = 2. for i = 0,5 do v(2,i) = 2. v(3,3) = 3. for i = 0,5 do v(3,i) = 3. v(4,4) = 4. for i = 0,5 do v(4,i) = 4. v(5,5) = 5. for i = 0,5 do v(5,i) = 5. print, 'testing str_field' if !version.os eq 'MacOS' then begin device, get_screen = ssize xsize = ssize(0)*.7 ysize = ssize(1)*.8 window, xsize = xsize, ysize = ysize end loadct, 12 print, 'Testing Vector Field Streamline approach' str_field1, u, v, xmargin=[5,5],ymargin=[0,0] stop, 'Press .con to continue..' print, 'Testing Level Sets of Stream Function approach' str_field2,u,v,x,y,xmargin=[5,5],ymargin=[5,5],/nolabels end ;********************************************************