pro boiler_plate, item, arg1 ;--------+---------+---------+---------+---------+---------+---------+-+ ; item = 1 Set plot device ; 2 Number of data files ; 3 Single input file ; 4 Multiple input files ; 5 Plot symbols ; 6 Special symbols ; 7 Output file name ; 8 Plot attributes ; 9 Close plot ; 10 Change device and replot ;--------+---------+---------+---------+---------+---------+---------+-+ ein = '($," Enter item number")' pcr = '($,/," Press to continue")' dataf = ' ' ans = ' ' if(item eq 2) then goto, j20 if(item eq 3) then goto, j30 if(item eq 4) then goto, j40 if(item eq 5) then goto, j50 if(item eq 6) then goto, j60 if(item eq 7) then goto, j70 if(item eq 8) then goto, j80 if(item eq 9) then goto, j90 if(item eq 10) then goto, j100 ;--------+ DEVICE TYPE j10: iws = 1 print, format='(/," Device type:")' print, format='(" 1) X-Window")' print, format='(" 2) Tektronix")' print, format='(" 3) PostScript File")' print, format=ein read, iws ; if (iws eq 1) then set_plot, 'sun' if (iws eq 1) then set_plot, 'X' if (iws eq 2) then set_plot, 'tek' if (iws eq 3) then begin set_plot, 'ps' print, format='($,/," Enter name of PostScript file")' read, dataf device, filename=dataf ; print, format='($,/," Insider trick? [N]")' ; read, ans ; if (strupcase(ans) eq 'Y') then $ ; device,output='85 35 {dup mul exch dup mul add 1.0 exch sub} setscreen' print, format='($,/," Encapsulated? [N]")' read, ans if (strupcase(ans) eq 'Y') then device, /encapsulated device, /inches, xoffset=0.0, yoffset=1.25 endif arg1 = iws return ;--------+ Number of input files j20: ndf = 1 print, format='($,/," Enter number of input data files")' read, ndf arg1 = ndf return ;--------+ Single input file j30: print, format='($,/," Enter name of input file")' read, dataf openr, unit, dataf, /get_lun arg1 = unit return ;--------+ Multiple input files j40: print, arg1, format='($,/," Enter name of input file no. ",i2)' read, dataf openr, unit, dataf, /get_lun arg1 = unit return ;--------+ Plot symbol or line type j50: ign = 1 print, format='(/," Plot symbol type:")' print, format='(" 1) Plus sign ",T30," 2) X")' print, format='(" 3) Asterisk ",T30," 4) Circle")' print, format='(" 5) Square ",T30," 6) Triangle")' print, format='(" 7) Diamond ",T30," 8) Hexagon")' print, format='(" 9) USD Triangle")' print, format=ein read, ign if (ign ge 4) then begin print, format='($,/," Fill symbol? [N]")' read, ans if (strupcase(ans) eq 'Y') then ign = -ign endif arg1 = ign return ;--------+ Special symbols j60: ign = arg1 if (abs(ign) eq 1) then begin ;Plus x = [-1,1,0,0,0] & y = [0,0,0,1,-1] endif if (abs(ign) eq 2) then begin ;X x = [-1,1,0,-1,1] & y = [-1,1,0,1,-1] endif if (abs(ign) eq 3) then begin ;Asterisk a = sin(!pi/4.0) x = [-1,1,0,0,0,0,-a,a,0,-a,a] & y = [0,0,0,1,-1,0,-a,a,0,a,-a] endif if (abs(ign) eq 4) then begin ;Circle a = findgen(17)*(!pi*2.0/16.0) x = cos(a) & y = sin(a) endif if (abs(ign) eq 5) then begin ;Square x = [-1,1,1,-1,-1] & y = [-1,-1,1,1,-1] endif if (abs(ign) eq 6) then begin ;Triangle a = 2.0*tan(!pi/6.0) & b = 2.0/3.0 x = [0,-a,a,0] & y = [2.0*b,-b,-b,2.0*b] endif if (abs(ign) eq 7) then begin ;Diamond a = sqrt(2.0) x = [0,-a,0,a,0] & y = [a,0,-a,0,a] endif if (abs(ign) eq 8) then begin ;Hexagon a = tan(!pi/6.0) x = [2.0*a,a,-a,-2.0*a,-a,a,2.0*a] & y = [0,1,1,0,-1,-1,0] endif if (abs(ign) eq 9) then begin ;USD Triangle a = 2.0*tan(!pi/6.0) & b = 2.0/3.0 x = [-a,0,a,-a] & y = [b,-2.0*b,b,b] endif usersym, x, y, fill=(ign lt 0) return ;--------+ Output file j70: print, format='($,/," Enter name of output file")' read, dataf arg1 = dataf return ;--------+ Change plot attributes j80: iws = arg1 ; if (iws le 2) then begin ; print, format='(/," Do you want to:")' ; print, format='(" 1) Change plot attributes and replot")' ; print, format='(" 2) Return to main menu [Default]")' ; print, format=ein ; read, ans ; if (ans eq '1') then begin ; plot_attrib, xc, yc, xshift, yshift, xtn, ytn ; if (!d.name eq 'SUN') then wdelete ; if (!d.name eq 'TEK') then erase ; goto, j400 ; endif ; if (!d.name eq 'SUN') then wdelete ; if (!d.name eq 'TEK') then erase ; endif return ;--------+ Close plot j90: iws = arg1 if (iws le 2) then begin print, format=pcr read, ans if (!d.name eq 'TEK') then erase if (!d.name eq 'X') then wdelete endif return ;--------+ Change device and replot j100: print, format='($,/," Change device and replot? [N]")' read, ans if (strupcase(ans) eq 'Y') then goto, j10 arg1 = 0 return end pro cntr_read, y, z, u ;--------+---------+---------+---------+---------+---------+---------+-+ ;READ CONTOUR FIELD ;--------+---------+---------+---------+---------+---------+---------+-+ common parm1, ids, title datin = 'um_5.wav' title = ' ' ids = 1 ny = 1 nz = 1 print, format='($,/," Enter contour input data file name")' read, datin openr, un1, datin, /get_lun readf, un1, ids readf, un1, title, format='(1a72)' readf, un1, ny, nz ;DECLARE DATA TYPE AND DIMENSION dummy = ' ' yd = fltarr(ny) y = fltarr(ny,nz) z = fltarr(nz,ny) u = fltarr(nz,ny) ;READ DATA readf, un1, yd readf, un1, z readf, un1, u free_lun, un1 for k=0,nz-1 do y(*,k) = yd z = transpose(z) u = transpose(u) return end pro plot_size, x_wnd,y_wnd,x_reg,y_reg ;--------+---------+---------+---------+---------+---------+---------+-+ ; ;--------+---------+---------+---------+---------+---------+---------+-+ if (!p.charsize eq 0) then !p.charsize=1.0 iend = 1 if (!d.name eq 'PS') then iend=5 for i=1,iend do begin xlm = !x.margin(0)*!d.x_ch_size*!p.charsize/(!d.x_px_cm*2.54) xrm = !x.margin(1)*!d.x_ch_size*!p.charsize/(!d.x_px_cm*2.54) ybm = !y.margin(0)*!d.y_ch_size*!p.charsize/(!d.y_px_cm*2.54) ytm = !y.margin(1)*!d.y_ch_size*!p.charsize/(!d.y_px_cm*2.54) x_reg = xlm+x_wnd+xrm y_reg = ybm+y_wnd+ytm ; print, !d.x_ch_size, !d.y_ch_size ; print, x_reg, y_reg if (!d.name eq 'X' or !d.name eq 'SUN') then begin x_reg = fix(x_reg*(!d.x_px_cm*2.54)) y_reg = fix(y_reg*(!d.y_px_cm*2.54)) window, xsize=x_reg, ysize=y_reg endif if (!d.name eq 'PS') then begin device, xsize=x_reg, ysize=y_reg, /inches endif endfor return end pro vect_read, y, z, v, w ;--------+---------+---------+---------+---------+---------+---------+-+ ;READ VECTOR FIELD ;--------+---------+---------+---------+---------+---------+---------+-+ common parm1, ids, title datin = 'vr_5.wav' title = ' ' ids = 1 ny = 1 nz = 1 print, format='($,/," Enter vector input data file name")' read, datin openr, un1, datin, /get_lun readf, un1, ids readf, un1, title, format='(1a72)' readf, un1, ny, nz ;DECLARE DATA TYPE AND DIMENSION dummy = ' ' y = fltarr(ny) zd = fltarr(nz,ny) v = fltarr(nz,ny) w = fltarr(nz,ny) ;READ DATA readf, un1, y readf, un1, zd readf, un1, v readf, un1, w free_lun, un1 z = zd(*,0) v = transpose(v) w = transpose(w) return end pro vect, u, v, x, y, sfact ;--------+---------+---------+---------+---------+---------+---------+-+ ; Plot vector field ;--------+---------+---------+---------+---------+---------+---------+-+ s = size(u) t = size(v) mag = sqrt(u^2+v^2) ;magnitude. maxmag = max(mag)/sfact x0 = min(x) ;get scaling x1 = max(x) y0 = min(y) y1 = max(y) sina = (x1-x0)/s(1)/maxmag*u ;sin & cosine components. cosa = (y1-y0)/s(2)/maxmag*v ; psymsav=!p.psym !p.psym=3 ;plot points ; plot,[x0-1.,x1+1.],[y1+1.,y0-1.] ;make grid. ; r = .3 ;len of arrow head st = r * 0.38268 ;sin 22.5 degs * length of head ct = r * 0.92387 ; for iy=0,s(2)-1 do for ix=0,s(1)-1 do begin if mag(ix,iy) gt 0.001 then begin x0 = x(ix) ;get coords of start & end dx = sina(ix,iy) x1 = x0 + dx y0 = y(iy) dy = cosa(ix,iy) 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)];, /data endif endfor ; oplot,[x0-1.,x1+1.],[y1+1.,y0-1.] ;make grid. !p.psym=psymsav return end pro wall, y, z ;--------+---------+---------+---------+---------+---------+---------+-+ ;DEFINE DUCT WALLS ;--------+---------+---------+---------+---------+---------+---------+-+ common parm1, ids, title ;DATA STATION rad = 10.214 a = [1.0,1.0,1.1734,1.3729,1.5464,1.5464] b = [1.0,1.0,0.8462,0.6693,0.5154,0.5154] eta = [2.0,2.0,3.4528,4.8791,10.000,10.000] a = a*rad b = b*rad npts = 101 ang = fltarr(npts) y = fltarr(npts) z = fltarr(npts) h = !pi/(2.0*float(npts-1)) ang(0) = 0.0 y(0) = a(ids)/rad z(0) = 0.0 for i=1,npts-2 do begin ang(i) = ang(i-1)+h cos1 = cos(ang(i)) sin1 = sin(ang(i)) r = (1.0/((cos1/a(ids))^eta(ids) + $ (sin1/b(ids))^eta(ids)))^(1.0/eta(ids)) y(i) = r*cos1/rad z(i) = r*sin1/rad endfor y(npts-1) = 0.0 z(npts-1) = b(ids)/rad ; yref = a/2.0 ; zref = b*(1.0-(yref/a)^eta)^(1.0/eta)+1.0 return end ;--------+---------+---------+---------+---------+---------+---------+-+ ; *** MAIN PROGRAM *** ; ;--------+---------+---------+---------+---------+---------+---------+-+ common parm1, ids, title ans = ' ' dummy = ' ' ein = '($," Enter item number")' a = [1.0,1.0,1.1734,1.3729,1.5464,1.5464] b = [1.0,1.0,0.8462,0.6693,0.5154,0.5154] yrn = [1.1,1.1,1.0,0.8,0.6,0.6] skey = [' ',' ', $ '!5V!Dr!N/U!Db!N=0.32','!5V!Dr!N/U!Db!N=0.32', $ '!5V!Dr!N/U!Db!N=0.12','!5V!Dr!N/U!Db!N=0.12'] rad = 10.214 ;--------+ DEVICE TYPE boiler_plate, 1, iws ;--------+ PLOT TYPE ipt = 1 print, format='(/," Plot type:")' print, format='(" 1) Contour plot")' print, format='(" 2) Vector plot")' print, format=ein read, ipt ;--------+ READ DATA if (ipt eq 1) then cntr_read, yax, zax, u if (ipt eq 2) then vect_read, yvr, zvr, v, w ;--------+ CONSTRUCT WALL ids = ids-1 wall, yw, zw xcl = [0.0] ycl = [0.0] a1 = a(ids) b1 = b(ids) ac = a1/2.0 ;--------+ PLOT! !x.thick=2.0 ;double axis thickness !y.thick=2.0 !p.charsize = 1.4 ; title = '!5' + strtrim(title) title = '!5' xrng = [-0.1,1.6] xspn = xrng(1)-xrng(0) xwid = 5.0 yrng = [-yrn(ids),yrn(ids)] yspn = yrng(1)-yrng(0) yhgt = xwid*(yspn/xspn) zrng = [0,10] j500: ;--------+ PLOT CONTOURS if (ipt eq 1) then begin cmin = min(u) cmax = max(u) print, format='(/," Min value = ",f10.5,5x,/," Max value = ",f10.5)', $ cmin, cmax print, format='($,/," Enter min contour level")' read, xmin print, format='($,/," Enter max contour level")' read, xmax print, format='($,/," Enter del contour level")' read, xdel nl = fix((xmax-xmin)/xdel) + 1 cl = xmin+findgen(nl)*xdel ;--------+ PLOT 2-D CONTOURS iplt = 'N' print, format='($,/," Plot 2-D contours? [N]")' read, iplt if (strupcase(iplt) eq 'Y') then begin plot_size, xwid,yhgt,dum1,dum2 plot, yw, zw, thick=3.0, $ title=title, $ xtitle='!5y/R', ytitle='!5z/R', $ xrange=xrng, yrange=yrng, xstyle=1, ystyle=1 oplot, yw, -zw, thick=3.0 !p.thick=1.0 oplot, [-0.03,0.03], [0.0,0.0] oplot, [0.0,0.0], [-0.03,0.03] oplot, [0.06,ac-0.06], [0.0,0.0] oplot, [ac-0.03,ac+0.03], [0.0,0.0] oplot, [ac+0.06,a1-0.03], [0.0,0.0] oplot, [0.0,0.0], [0.06,b1-0.03] oplot, [0.0,0.0], [-0.06,-b1+0.03] if (ids eq 2) then begin ; Peak vorticity ; boiler_plate, 6, -4 ; oplot, [1.04445,1.04445], [-0.541742,0.541742], psym=8 xyouts, 0.2,0.3,'(IMAGED)', size=1.4 endif if (ids eq 3) then begin ; boiler_plate, 6, -4 ; oplot, [1.29313,1.29313], [-0.419459,0.419459], psym=8 xyouts, 0.2,0.2,'(IMAGED)', size=1.4 endif ; if (ids eq 4) then begin ; boiler_plate, 6, -4 ; oplot, [1.44234,1.44234], [-0.0958483,0.0958483], psym=8 ; oplot, [0.0,1.5456], [0.1,0.1] ; oplot, [0.0,1.5456], [0.2,0.2] ; oplot, [0.0,1.5456], [0.3,0.3] ; oplot, [0.0,1.54], [0.4,0.4] ; oplot, [0.0,1.5456], [-0.1,-0.1] ; oplot, [0.0,1.5456], [-0.2,-0.2] ; oplot, [0.0,1.5456], [-0.3,-0.3] ; oplot, [0.0,1.54], [-0.4,-0.4] ; endif ; if (ids eq 5) then begin ; boiler_plate, 6, -4 ; oplot, [1.34287,1.34287], [-0.1,0.1], psym=8 ; oplot, [0.0,1.5456], [0.1,0.1] ; oplot, [0.0,1.5456], [0.2,0.2] ; oplot, [0.0,1.5456], [0.3,0.3] ; oplot, [0.0,1.54], [0.4,0.4] ; oplot, [0.0,1.5456], [-0.1,-0.1] ; oplot, [0.0,1.5456], [-0.2,-0.2] ; oplot, [0.0,1.5456], [-0.3,-0.3] ; oplot, [0.0,1.54], [-0.4,-0.4] ; endif ; if (ids eq 4 or ids eq 5) then begin ; y1 = yax(*,0) ; is = where(y1 ge 1.0, count) ; j1 = is(0)-1 ; is = size(yax) ; j2 = is(1) ; ns = j2-j1 ; yrdg = fltarr(ns) ; zrdg = yrdg ; for j=0,ns-1 do begin ; jdum = j+j1-1 ; yrdg(j) = yax(jdum,0) ; udum = u(jdum,*) ; zdum = zax(jdum,*) ; umax = max(udum) ; is = where(udum eq umax) ; zrdg(j) = zdum(is) ; endfor ; oplot, yrdg, zrdg, linestyle=2 ; endif !p.thick=2.0 contour, u, yax, zax, $ levels=cl, c_linestyle=2*(cl lt 0.0), $ xrange=xrng, yrange=yrng, xstyle=1, ystyle=1, $ /follow, /spline, /noerase print, format='($,/," Change contour levels and replot [N]")' read, ans if (strupcase(ans) eq 'Y') then begin erase goto, j500 endif endif ;--------+ PLOT 3-D SHADED SURFACE WITH WIRE MESH AND CONTOURS iplt = 'N' print, format='($,/," Plot shaded surface? [N]")' read, iplt if (strupcase(iplt) eq 'Y') then begin print, format='($,/," Enter z-axis maximum")' read, zr2 zrng = [0,zr2] ; s1 = (u ne 0) ; uref = u ; is = where(u eq 0,count) ; uref(is) = 1.0 ; u = uref plot_size, xwid,yhgt,dum1,dum2 surface, u, yax, zax, $ xrange = xrng, yrange = yrng, zrange = zrng, $ xstyle=5, ystyle=5, zstyle=5, /save, /nodata ; shade_surf, u, yax, zax, /t3d, $ ; xrange = xrng, yrange = yrng, zrange = zrng, $ ; xstyle=5, ystyle=5, zstyle=5, /noerase ; shade_surf, uref, yax, zax, /t3d, $ ; xrange = xrng, yrange = yrng, zrange = zrng, $ ; xstyle=5, ystyle=5, zstyle=5, /noerase oplot, yw, zw, thick=3.0, /t3d, zvalue=0, /noerase, /noclip oplot, yw,-zw, thick=3.0, /t3d, zvalue=0, /noerase, /noclip oplot, yw, zw, thick=3.0, /t3d, zvalue=1, /noerase, /noclip oplot, yw,-zw, thick=3.0, /t3d, zvalue=1, /noerase, /noclip oplot,[-0.03,0.03], [0.0,0.0],thick=1.0,/t3d,zvalue=1,/noerase,/noclip oplot,[0.0,0.0], [-0.03,0.03],thick=1.0,/t3d,zvalue=1,/noerase,/noclip oplot,[0.06,ac-0.06],[0.0,0.0],thick=1.0,/t3d,zvalue=1,/noerase,/noclip oplot,[ac-.03,ac+.03], [.0,.0],thick=1.0,/t3d,zvalue=1,/noerase,/noclip oplot,[ac+.06,a1-0.03],[.0,.0],thick=1.0,/t3d,zvalue=1,/noerase,/noclip oplot,[0.0,0.0],[0.06,b1-0.03],thick=1.0,/t3d,zvalue=1,/noerase,/noclip oplot,[.0,.0],[-0.06,-b1+0.03],thick=1.0,/t3d,zvalue=1,/noerase,/noclip surface, u, yax, zax, /t3d, skirt=0, $ xrange = xrng, yrange = yrng, zrange = zrng, $ ztitle = '!5', $ xtitle = '!5y/R', ytitle = '!5z/R', $ xstyle=1, ystyle=1, zstyle=1, /noerase ; surface, uref, yax, zax, /t3d, $ ; xrange = xrng, yrange = yrng, zrange = zrng, $ ; ztitle = '!5', $ ; xtitle = '!5y/R', ytitle = '!5z/R', $ ; xstyle=1, ystyle=1, zstyle=1, /noerase contour, /follow, u, yax, zax, /t3d, zvalue = 1, $ levels=cl, c_linestyle=2*(cl lt 0.0), $ xrange = xrng, yrange = yrng, $ ; xticks=6, xtickv=[-3,-2,-1,0,1,2,3], xminor=2, $ ; yticks=5, yminor=2, $ xstyle=1, ystyle=1, zstyle=1, /noerase ; cl = [1,2] ; contour, /follow, s1, yax, zax, /t3d, zvalue = 0.1, $ ; levels=cl, c_linestyle=2*(cl lt 0.0), $ ; xrange = xrng, yrange = yrng, thick=3.0, $ ; xstyle=5, ystyle=5, zstyle=5, /noerase endif endif ;--------+ PLOT VECTORS if (ipt eq 2) then begin sfact = 1.0 j600: plot_size, xwid,yhgt,dum1,dum2 print, sfact, format='($,/," Vector scaling = ",f3.1," Change [N]")' read, ans if (strupcase(ans) eq 'Y') then begin print, format='($,/," Enter vector scaling")' read, sfact endif plot, yw, zw, thick=3.0, $ title=title, $ xtitle='!5y/R', ytitle='!5z/R', $ xrange=xrng, yrange=yrng, xstyle=1, ystyle=1 oplot, yw, -zw, thick=3.0 !p.thick=1.0 oplot, [-0.03,0.03], [0.0,0.0] oplot, [0.0,0.0], [-0.03,0.03] xyouts, 1.2/4.02125, 0.57/4.02125, skey(ids) if (ids eq 4) then begin boiler_plate, 6, -4 oplot, [1.44234,1.44234], [-0.0958483,0.0958483], psym=8 endif if (ids eq 5) then begin boiler_plate, 6, -4 oplot, [1.34287,1.34287], [-0.1,0.1], psym=8 endif !p.thick=2.0 vect, v,w,yvr,zvr,sfact print, format='($,/," Change vector scaling and replot [N]")' read, ans if (strupcase(ans) eq 'Y') then begin erase goto, j600 endif endif ;--------+ Close plot if (iws le 2) then boiler_plate, 9, iws ;--------+ Change device and replot? boiler_plate, 10, iws if (iws ne 0) then goto, j500 if (!d.name eq 'PS') then device, /close_file stop, ' ' end