pro plot_vectmag, bx, by, bz, imgi, xres = xres, yres = yres, level = level, $ text = text, btcut = btcut, device = device, xrange=xrange, yrange=yrange, $ image = image, jpgfilename = jpgfilename, psfilename = psfilename, tcol = tcol, $ notrans = notrans ;+ ;PURPOSE: ; plot the vector magnetic field either on the screen, or in JPG, or in PostScript. ;INPUT: ; bx: x-component of transverse field. ; by: y-component of transverse field. ; bz: z-component, longitudinal field. ;OUTPUT: ; no explicit output; create screen/JPG/PS display. ;OPTIONAL INPUT: ; xres: image scale in X dimension, default = 315./512 (arcsec/pixel) ; yres: image scale in Y dimension, default = 315./512 (arcsec/pixel) ; level: the contour levels for Bz; default = [-0.8, -0.6, -0.4, -0.2, -0.1, $ ; 0.1, 0.2, 0.4, 0.6, 0.8]*max(abs(Bz)) ; text: figure caption you want to write on the image; if not set, keyboard input. ; btcut: the cutoff value of transverse field (the field below this value is not plotted); ; default = mean(Bt) ; device: specify output device; ; device = 0: screen output (default); user is then prompted to choose PS or JPG or none. ; device = 1: create JPG file without screen display; user prompted to choose PS. ; device = 2: create PS file. ; xrange: the X-range of sub-frame chosen to be displayed; xrange = [x1, x2] ; yrange: the Y-range of sub-frame chosen to be displayed; yrange = [y1, y2] ; image: specify which image (Bz or imgi) to use as background image ; image = 0: Bz as background ; image = 1: imgi as background ; jpgfilename: the output JPG filename; default = 'vectmag.jpg' ; psfilename: the output PS filename; default = 'vectmag.ps' ; tcol: color of the text. ;OPTIONAL OUTPUT: ; ;KEYWORDS: ; notrans: if set, no transverse vectors should be drawn, e.g., for quiet regions. ;HISTORY: ; 1998/06/11: first version at BBSO. QiuJ ; 2001/02/24: modifications made for Tom's display. QiuJ. ; 2001/06/07: modified a bit more to choose contour levels based on sub-frame image, ; and better adjust the vector plot grids, ; and add /notrans keyword for quiet region plot. QiuJ. ; 2001/09/05: change the output GIF into JPG file to be compatible with IDL_5.4 ;- szx=size(bx) szy=size(by) szz=size(bz) if szx(0) ne 2 or szy(0) ne 2 or szz(0) ne 2 then begin print,'Data must be 2D array.' goto, bail endif for i=0, 2 do begin if szx(i) ne szy(i) or szx(i) ne szz(i) or szy(i) ne szz(i) then begin print,'Image dimensions unmatch.' goto, bail endif end nx0=szz(1) ny0=szz(2) ; image scale if n_elements(xres) eq 0 then xres = 315./512. ; arcsec/pixel if n_elements(yres) eq 0 then yres = 315./512. ; arcsec/pixel ;maxbx=max(bx) & minbx=min(bx) ;maxby=max(by) & minby=min(by) ;maxbz=max(bz) & minbz=min(bz) ;print, ' ', 'Max', 'Min' ;print, 'bx ', maxbx, minbx ;print, 'by ', maxby, minby ;print, 'bz ', maxbz, minbz if n_elements(text) eq 0 then begin date='' time='' numb='' print, '' read,'Input the date of the magnetogram: ', date read,'Input the time of the magnetogram: ', time read,'Input the AR # of the magnetogram: ', numb print, '' text = date + ' ' + time + ' NOAA ' + numb endif if n_elements(xrange) eq 0 or n_elements(yrange) eq 0 then begin yn='' print, '' read,'Do you want to extract sub-frame (y/n) ==> ', yn if yn ne 'y' and yn ne 'Y' then begin xrange=[0, nx0-1] yrange=[0, ny0-1] goto, draw endif else begin set_plot,'x' window, 1, xs = nx0, ys = ny0, r = 2 region, bz, xy, /norep, /nonum xrange=[xy(0,0), xy(0,1)] yrange=[xy(1,0), xy(1,1)] endelse endif draw: x1=xrange(0) & x2=xrange(1) & y1=yrange(0) & y2=yrange(1) print,'Sub-frame X-range: ('+string(x1, format = '(i3)') + ', ' + string(x2, format = '(i3)') +')' print,'Sub-frame y-range: ('+string(y1, format = '(i3)') + ', ' + string(y2, format = '(i3)') +')' print, '' nx=abs(x2-x1+1) & ny=abs(y2-y1+1) xct=indgen(nx)*xres yct=indgen(ny)*yres ; contour levels for Bz; when the calibration is done, one can change it to Gauss. if n_elements(level) eq 0 then level = [-0.8, -0.6, -0.4, -0.2, -0.1, 0.1, 0.2, 0.4, 0.6, 0.8] * max(abs(bz(x1:x2, y1:y2))) level = level(sort(level)) ; contour levels must be in order n_leveln = n_elements(where(level lt 0)) ; how many negative contours n_levelp = n_elements(where(level gt 0)) ; how many positive contours ; here set a transverse field cutoff value; if n_elements(btcut) eq 0 then begin bt = sqrt(bx(x1:x2, y1:y2)^2+by(x1:x2, y1:y2)^2) ;sigma = stdev(Bt) btcut = mean(bt) print, 'Transverse Field (Bt) Cutoff ==' + string(btcut, format = '(f5.2)') print, '' endif if keyword_set(notrans) then notr = 1 else notr = 0 ; ok make the color map: r = bindgen(256) b = bindgen(256) g = bindgen(256) ; red b(254) = 0 g(254) = 0 ; blue r(253) = 0 g(253) = 0 ; green r(252) = 0 b(252) = 100 ;; or if you like, yellow ;b(252) = 0 tvlct, r, g, b c_color = [replicate(253, n_leveln), replicate(254, n_levelp)] if n_elements(image) eq 0 then image = 0 if image then bac_img = bytscl(imgi0.2*max(imgi), top = 251) else $ bac_img = bytscl(bz<0.5*max(abs(bz))>(-0.5)*max(abs(bz)), top = 251) if n_elements(tcol) eq 0 then if image then tcol = 0 else tcol = 255 ; what device for display? if n_elements(device) eq 0 then read, 'Screen[0] or JPG[1] or PS[2] ==> ', device print, '' dev=device psyn='' gfyn='' if (dev eq 0) then begin ; screen display set_plot, 'x' zmx=640/nx zmy=640/ny zm=(zmx gt zmy)*zmy+(zmx le zmy)*zmx x_siz = zm*nx & y_siz=zm*ny window, 0, xsize = x_siz+60, ysize = y_siz+60, r = 2 tv, rebin(bac_img(x1:x2, y1:y2), x_siz, y_siz), 40, 40 contour, bz(x1:x2, y1:y2), xct, yct, pos=[40,40,x_siz+39, y_siz+39], xst=1, yst=1,$ /dev, level = level, title=title, c_thick = 2, c_color = c_color, xtitle = '!17arcsecond',$ ytitle = '!17arcsecond', /noerase contour, smooth(bz(x1:x2, y1:y2),4, /edge), xct, yct, pos=[40,40,x_siz+39, y_siz+39], xst=1, yst=1,$ /dev, level=[0],/noerase, c_color = 0, c_thick = 3 if not notr then vect1, bx(x1:x2, y1:y2), by(x1:x2, y1:y2), 1, nx/4, ny/4, rlen = 0.4, color = 252, thick = 2, $ pos=[40,40,x_siz+39, y_siz+39],xst=5, yst=5, missing=btcut, /dev, /noer, length = 1.2 xyouts, 50, y_siz + 10, '!5 BBSO Vector Magnetogram', size = 2, /dev, color = tcol xyouts, 50, y_siz - 10, '!5 '+text, size = 2, /dev, color = tcol read, 'Make JPG file(y/n) ==> ', gfyn read, 'Make PS file(y/n) ==> ', psyn print, '' endif ; screen display if (dev eq 1 or gfyn eq 'y' or gfyn eq 'Y') then begin ; make JPG file set_plot, 'z' zmx=640/nx zmy=640/ny zm=(zmx gt zmy)*zmy+(zmx le zmy)*zmx x_siz = zm*nx & y_siz=zm*ny device, set_resolution = [x_siz+60, y_siz+60] tv, rebin(bac_img(x1:x2, y1:y2), x_siz, y_siz), 40, 40 contour, bz(x1:x2, y1:y2), xct, yct, pos=[40,40,x_siz+39, y_siz+39], xst=1, yst=1,$ /dev, level = level, c_thick = 2, c_color = c_color, xtitle = '!17arcsecond', $ ytitle = '!17arcsecond', /noerase contour, smooth(bz(x1:x2, y1:y2),4, /edge), xct, yct, pos=[40,40,x_siz+39, y_siz+39], xst=1, yst=1,$ /dev, level=[0], /noerase, c_color = 0, c_thick = 3 if not notr then vect1, bx(x1:x2, y1:y2), by(x1:x2, y1:y2), 1, nx/4, ny/4, rlen = 0.4, color = 252, thick = 2, $ pos=[40,40,x_siz+39, y_siz+39],xst=5, yst=5, missing=btcut, /dev, /noer, length = 1.2 xyouts, 50, y_siz + 10, '!5 BBSO Vector Magnetogram', size = 2, /dev, color = tcol xyouts, 50, y_siz - 10, '!5 '+text, size = 2, /dev, color = tcol t = tvrd() ;if n_elements(giffilename) eq 0 then giffilename = 'vectmag.jpg' if n_elements(jpgfilename) eq 0 then jpgfilename = 'vectmag.jpg' ;print, 'Write GIF file into ' + giffilename print, 'Write JPEG file into ' + jpgfilename print, '' ;write_gif, giffilename, t, r, g, b write_image, jpgfilename, 'JPEG', t, r, g, b endif ; make JPG file if (dev eq 2 or psyn eq 'y' or psyn eq 'Y') then begin ; make PS file if n_elements(psfilename) eq 0 then psfilename = 'vectmag.ps' print, 'Write PS file into ' + psfilename print, '' if nx gt ny then pt=0 else pt=1 if n_elements(pssize) eq 0 then pssize = 16. if pt eq 0 then begin psxsize=pssize + 3. psysize=(psxsize - 3)/float(nx)*float(ny) + 3. endif else begin psysize=pssize + 3. psxsize=(psysize - 3)/float(ny)*float(nx) + 3. endelse set_plot,'ps' if pt eq 0 then device, filename=psfilename, xsize=psxsize, ysize=psysize, $ /landscape, /color, /bold, /times, /iso, font_index = 4, bits = 8 if pt eq 1 then device, filename=psfilename, xsize=psxsize, ysize=psysize, $ xoffset=(18-psxsize)/2., yoffset=(25.2-psysize)/2.,/portrait, /color, $ /bold, /times, /iso, font_index = 4, bits = 8 ;now draw the things tv, bac_img(x1:x2, y1:y2), 2., 2., xs = psxsize -3., ys = psysize -3., /cent contour, bz(x1:x2, y1:y2), xct, yct, pos=[2/psxsize, 2/psysize, (psxsize-1.)/psxsize, $ (psysize - 1.)/psysize], xst=1, yst=1,/normal, /noerase, $ level = level, c_color = c_color, thick = 4, xtitle = '!4arcsecond', ytit = '!4arcsecond', font = 0 contour, smooth(bz(x1:x2, y1:y2),4, /edge), xct, yct, pos=[2./psxsize, 2./psysize, (psxsize-1.)/psxsize,$ (psysize - 1.)/psysize], xst=5, $ yst=5,/normal, level = 0 ,c_color = 0, c_thick = 5, /noerase if not notr then vect1, bx(x1:x2, y1:y2), by(x1:x2, y1:y2), 1, nx/4, ny/4, rlen = 0.4, thick = 5, $ pos=[2./psxsize, 2./psysize, (psxsize-1.)/psxsize, (psysize - 1.)/psysize],xst = 5, yst = 5, $ /nor, /noer, missing=btcut, color = 252, length = 1.2 xyouts, 2.5/psxsize, (psysize - 2.)/psysize, '!4 BBSO Vector Magnetogram', color = tcol, font = 0, size = 1.5, /nor xyouts, 2.5/psxsize, (psysize - 3.)/psysize, '!4 '+text, font = 0, size = 1.5, /nor, color = tcol device,/close endif ; make PS file set_plot, 'x' ; get back to screen set. bail: end ;********************************************************************** PRO VECT1,UU,VV,ISARROW,SX,SY,X,Y, Missing = Missing, Length = length, $ Dots = dots, Color=color, Thick=thick, Rlen=rlen, T3d=t3d, $ _EXTRA = extra ; ;+ ; NAME: ; VELOVECT ; ; PURPOSE: ; Produce a two-dimensional velocity field plot. ; ; A directed arrow is drawn at each point showing the direction and ; magnitude of the field. ; ; CATEGORY: ; Plotting, two-dimensional. ; ; CALLING SEQUENCE: ; VELOVECT, 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 smaller ; 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. ; ; COLOR: The color index used for the plot. ; ; RLEN: The length of arrow head. If rlen=0.0, then no arrow. ; ; Note: All other keywords are passed directly to the PLOT procedure ; and may be used to set option such as TITLE, POSITION, ; NOERASE, etc. ; OUTPUTS: ; None. ; ; COMMON BLOCKS: ; None. ; ; SIDE EFFECTS: ; Plotting on the selected device is performed. System ; variables concerning plotting are changed. ; ; RESTRICTIONS: ; None. ; ; PROCEDURE: ; Straightforward. Unrecognized keywords are passed to the PLOT ; procedure. ; ; MODIFICATION HISTORY: ; DMS, RSI, Oct., 1983. ; For Sun, DMS, RSI, April, 1989. ; Added TITLE, Oct, 1990. ; Added POSITION, NOERASE, COLOR, Feb 91, RES. ; August, 1993. Vince Patrick, Adv. Visualization Lab, U. of Maryland, ; fixed errors in math. ; August, 1993. DMS, Added _EXTRA keyword inheritance. ;- ; on_error,2 ;Return to caller if an error occurs u=uu v=vv s = size(uu) t = size(vv) 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 ; ;PRINT,'PARAMETER NUMBERS:',N_PARAMS() if n_params() lt 3 then isarrow=1. if n_params() lt 4 then SX=S(1) if n_params() lt 5 then SY=S(2) if n_params(0) lt 6 then x = findgen(SX) 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 7 then y = findgen(SY) else $ if n_elements(y) ne s(2) then goto,badxy ; if n_elements(missing) le 0 then missing = 1.0e-30 if n_elements(length) le 0 then length = 1.0 if n_elements(rlen) le 0 then rlen=isarrow*0.3 if n_elements(t3d) le 0 then t3d=0 if n_elements(thick) le 0 then thick=1 IF (SX NE S(1)) OR (SY NE S(2)) THEN BEGIN U=CONGRID(UU,SX,SY) V=CONGRID(VV,SX,SY) ENDIF 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 gt missing) if keyword_set(dots) then bad = where(mag le missing, nbad) endif else begin good = lindgen(n_elements(mag)) endelse ugood = u(good) vgood = v(good) x0 = min(x) ;get scaling x1 = max(x) y0 = min(y) y1 = max(y) x_step=(x1-x0)/sx y_step=(y1-y0)/sy maxmag=max([max(ugood/x_step),max(vgood/y_step)]) sina = length * (ugood/maxmag) cosa = length * (vgood/maxmag) ; if n_elements(title) le 0 then title = '' ;-------------- plot to get axes --------------- if n_elements(color) eq 0 then color = !p.color x_b0=x0-x_step x_b1=x1+x_step y_b0=y0-y_step y_b1=y1+y_step if n_elements(position) eq 0 then begin plot,[x_b0,x_b1],[y_b1,y_b0],/nodata,/xst,/yst, $ color=color, t3d=t3d, _EXTRA = extra endif else begin plot,[x_b0,x_b1],[y_b1,y_b0],/nodata,/xst,/yst, $ color=color, t3d=t3d, _EXTRA = extra endelse ; ; if keyword_set(arrow) then r=0.0 else r=.3 ; r = .3 ;len of arrow head r = rlen angle = 22.5 * !dtor ;Angle of arrowhead st = r * sin(angle) ;sin 22.5 degs * length of head ct = r * cos(angle) ; for i=0L,n_elements(good)-1 do begin ;Each point x0 = x(good(i) mod sx) ;get coords of start & end dx = sina(i) x1 = x0 + dx y0 = y(good(i) / sx) dy = cosa(i) y1 = y0 + dy xd=x_step yd=y_step plots,[x0,x1,x1-(ct*dx/xd-st*dy/yd)*xd, $ x1,x1-(ct*dx/xd+st*dy/yd)*xd], $ [y0,y1,y1-(ct*dy/yd+st*dx/xd)*yd, $ y1,y1-(ct*dy/yd-st*dx/xd)*yd], $ color=color,t3d=t3d, thick=thick endfor if nbad gt 0 then $ ;Dots for missing? oplot, x(bad mod sx), y(bad / sx), psym=3, $ color=color, t3d=t3d, thick=thick end ;************************************************************************************* pro region,img,xy,zoom=zoom,norep=norep,onlyplot=onlyplot,$ noplot1=noplot1,noplot2=noplot2,nonum=nonum,line=line,$ contour=contour,clevel=clevel,level=level,slevel=slevel,$ cont1=cont1,cont2=cont2,lev1=lev1,lev2=lev2,col1=col1,col2=col2, $ noimg=noimg ; ;+ ;NAME: ; region ;FUNCTION: ; interactively click out and show in the screen the regions in ; the image. ;INPUTS: ; img=arr(*,*) input image. ;OUTPUTS: ; xy=long(2,2,*): the downleft and upreght corners of the clicked ; regions. ; xy(0,0,*)= x of downleft corner ; xy(1,0,*)= y of downleft corner ; xy(0,1,*)= x of upright corner ; xy(1,1,*)= y of upright corner ;KEYWORDS: ; zoom: image display zoom factor ; norep: if set, click only one region ; nonum: if set do not remark the number ; onlyplot: plot out the regions included in the xy input ; without clicking the new ones. ; noplot1: if set, not plot the regions already in the xy but only ; those new ones. ; noplot2: if set, not plot the new clicked ones but only the old ; ones. ; contour: if set, make contours of the image. ; noimg: if set, not display the image. ;OPTIONAL INPUT: ; line: if input, draw the lines(such as the slits) on the image. ; cont1,cont2: other images used to overlay as contours, such as ; SXT and HXT images. ; lev1,lev2: contour levels for cont1 and cont2. ; clevel: contour level, if contour is to be made. ; slevel: the level for contouring the spot. ; col1,col2: the colors for cont1 and cont2. ; level: optional level used in drawline1.pro to draw the lines. ;OTHERS: ; this procedure will be used to make lightcurves. ;95/12/5 fixed Qiu J. Arcetri ;minor changes in 96. ;- ; if n_elements(zoom) eq 0 then zm=1 else zm=zoom if keyword_set(noimg) then s_tv=0 else s_tv=1 if s_tv eq 1 then sizm=size(img) ;tvscl,rebin(img,sizm(1)*zm,sizm(2)*zm) if keyword_set(contour) then drawline1,img,line,zoom=zm,level=level,$ slevel=slevel else if n_elements(line) ne 0 then drawline,img,line,zoom=zm $ ,color=220 else if s_tv eq 1 then tvscl,rebin(img,sizm(1)*zm,sizm(2)*zm) if n_elements(cont1) ne 0 then contour,cont1,pos=[0,0,sizm(1)*zm-1,$ sizm(2)*zm-1],/dev,/noerase,xst=1,yst=1,lev=lev1*max(cont1),c_color=col1 if n_elements(cont2) ne 0 then contour,cont2,pos=[0,0,sizm(1)*zm-1,$ sizm(2)*zm-1],/dev,/noerase,xst=1,yst=1,lev=lev2*max(cont2),c_color=col2 i=0 if n_elements(xy) ne 0 then if max(xy) gt 0 then begin sizxy=size(xy) if sizxy(0) eq 2 then begin sizxy=[sizxy(0)+1,sizxy(1),sizxy(2),1,sizxy(3),sizxy(4)] xy=reform(xy,2,2,1) endif nn=sizxy(sizxy(0)) if not keyword_set(noplot1) then begin for i=0,nn-1 do begin x1=xy(0,0,i)*zm y1=xy(1,0,i)*zm x2=xy(0,1,i)*zm y2=xy(1,1,i)*zm plots,[x1,x2],[y1,y1],/dev,color=255 plots,[x2,x2],[y1,y2],/dev,color=255 plots,[x2,x1],[y2,y2],/dev,color=255 plots,[x1,x1],[y2,y1],/dev,color=255 if not keyword_set(nonum) then begin ind=numname('#',i) xyouts,x2,y2,ind,/dev,size=1.5 endif endfor endif i=nn endif more: yn='y' if keyword_set(onlyplot) then goto,bail print, '' print,'Use your mouse to choose the sub-region:' print, '' print, 'Click to locate lower-left corner ==> ' print, '' cursor,x1,y1,/dev,/down print, 'Lower-left corner at: (' + string(x1, format = '(i3)')+ ', ' + string(y1, format = '(i3)')+')' print, '' print, 'Click to locate upper-right corner ==> ' print,'' cursor,x2,y2,/dev,/down print, 'Upper-right corner at: (' + string(x2, format = '(i3)')+ ', ' + string(y2, format = '(i3)')+')' print, '' if not keyword_set(noplot2) then begin plots,[x1,x2],[y1,y1],/dev plots,[x2,x2],[y1,y2],/dev plots,[x2,x1],[y2,y2],/dev plots,[x1,x1],[y2,y1],/dev if not keyword_set(nonum) then begin ind=numname('#',i) xyouts,x2,y2,ind,/dev,size=1.5 endif endif x1=x1/zm & y1=y1/zm & x2=x2/zm & y2=y2/zm if i eq 0 then xy=[[x1,y1],[x2,y2]] else xy=[[[xy]],[[x1,y1],[x2,y2]]] if not keyword_set(norep) then begin print,'more points?' read,yn if string(yn) eq 'n' then goto,bail i=i+1 goto,more endif bail: end