plot_psd.pro

pic 模拟程序！面向对象

    !p.font=-1
    !p.charsize=2.0
    device,decomposed=0
  endif

endif


; *****************************************************************
; Loop for rendering 2-D color contour plot on screen and to a file
; *****************************************************************

if (doColorCon eq 1) then begin

; Get a new window
  window_number = !d.window + 1
  print, ' '
  print, 'Color contour plot will appear in window ', window_number
  window, window_number
  ccontour_i = 0
  ccontour_jump:

; Load in the STD GAMMA-II color table
  loadct, 5

; Specify a font that looks great for printing (crappy on screen),
;   or else one that looks OK on the screen (also OK for printing).
  if (!d.name eq 'PS') then begin
    !p.font=1
    !p.charsize=1.6
    !p.charthick=1.5
    print, ' '
    print, 'Writing the 2D color contour plot to file ' + ps_ccon_file
  endif else begin
    !p.font=-1
    !p.charsize=2.0
    device,decomposed=0
    print, ' '
    print, 'Rendering the 2D color contour plot to the screen...'
  endelse

if (log_plot_flag eq 0) then begin
  contour, ez, z, r, thick=2, /xs, /ys, nlevels=100, xtitle=x_label, ytitle=y_label, /fill
endif else begin
  contour, ez_log, z, r, thick=2, /xs, /ys, nlevels=256, xtitle=x_label, ytitle=y_label, MIN_VALUE=28.0, MAX_VALUE=52.0, /fill
endelse

; rendering to screen is done.  here we render to ps file
  if (ccontour_i eq 0) then begin
    ccontour_i = 1
    set_plot, 'PS', /copy
    device, filename = ps_ccon_file, /color, bits_per_pixel=8
    goto, ccontour_jump  
  endif

; here rendering to file is done.  set back to screen defaults
  if (ccontour_i eq 1) then begin
    device, /close
    set_plot, 'X'
    !p.font=-1
    !p.charsize=2.0
    device,decomposed=0
  endif

endif


; *****************************************************************
; Loop for rendering color surface plot on screen and to a file
; *****************************************************************

if (doSurface eq 1) then begin

; Get a new window
  window_number = !d.window + 1
  print, ' '
  print, 'Color surface plot will appear in window ', window_number
  window, window_number
  surface_i = 0
  surface_jump:

; Load in the STD GAMMA-II color table
  loadct, 5

; Specify a font that looks great for printing (crappy on screen),
;   or else one that looks OK on the screen (also OK for printing).
  if (!d.name eq 'PS') then begin
    !p.font=1
    !p.charsize=1.6
    !p.charthick=1.5
    print, ' '
    print, 'Writing the color surface plot to file ' + ps_surf_file
  endif else begin
    !p.font=-1
    !p.charsize=2.0 
    device,decomposed=0
    print, ' '
    print, 'Rendering the color surface plot to the screen...'
  endelse

  xAng = 50;
  zAng = 30;

if (log_plot_flag eq 0) then begin
  surface, ez, z, r, thick=2, /xs, /ys, ax=xAng, az=zAng, shades=ez_color, $
           xtitle=x_label, ytitle=y_label, ztitle=z_label
endif else begin
  surface, ez, z, r, thick=2, /xs, /ys, ax=xAng, az=zAng, shades=ez_color, $
           xtitle=x_label, ytitle=y_label, ztitle=z_label, /zlog
endelse

; rendering to screen is done.  here we render to ps file
  if (surface_i eq 0) then begin
    surface_i = 1
    set_plot, 'PS', /copy
    device, filename = ps_surf_file, /color, bits_per_pixel=8
    goto, surface_jump  
  endif

; here rendering to file is done.  set back to screen defaults
  if (surface_i eq 1) then begin
    device, /close
    set_plot, 'X'
    !p.font=-1
    !p.charsize=2.0
    device,decomposed=0
  endif

endif


; *****************************************************************
; Repeat simple plot loop for generating fancy 3-tiered plot.
; *****************************************************************

if (doShow3D eq 1) then begin

; Get a new window
  window_number = !d.window + 1
  print, ' '
  print, 'Three-tiered plot will appear in window ', window_number
  window, window_number
  three_tiered_i = 0
  three_tiered_jump:

; Load in the STD GAMMA-II color table
  loadct, 5

; Specify a font that looks great for printing (crappy on screen),
;   or else one that looks OK on the screen (also OK for printing).
  if (!d.name eq 'PS') then begin
    !p.font=1
    !p.charsize=1.6
    !p.charthick=1.5
    print, ' '
    print, 'Writing 3-tiered surface/contour plot to file ' + ps_tier_file
  endif else begin
    !p.font=-1
    !p.charsize=2.0
    device,decomposed=0
    print, ' '
    print, 'Rendering the 3-tiered surface/contour plot to the screen...'
  endelse

  surface_struct = { shades:ez_color, xtitle:x_label, ytitle:y_label, ztitle:z_label, $
                     xstyle:1, ystyle:1, zstyle:1}
  contour_struct = { nlevels:6, xstyle:1, ystyle:1, zstyle:1}
  show3, ez, z, r, e_surface=surface_struct, e_contour=contour_struct

; rendering to screen is done.  here we render to ps file
  if (three_tiered_i eq 0) then begin
    three_tiered_i = 1
    set_plot, 'PS', /copy
    device, filename = ps_tier_file, /color, bits_per_pixel=8
    goto, three_tiered_jump  
  endif

; here rendering to file is done.  set back to screen defaults
  if (three_tiered_i eq 1) then begin
    device, /close
    set_plot, 'X'
    !p.font=-1
    !p.charsize=2.0
    device,decomposed=0
  endif

endif

; *****************************************************************
; Repeat simple plot loop for generating 2-D lineout plots.
; *****************************************************************

if (doLineout eq 1) then begin
;
; First, grab the desired column of data, the z array corresponds to
; the first coordinate and the r array to the second one, i.e. to
; rows and columns.
; Specify which column of grid points to use with the z_factor.
; Set the z_factor (between 0 and 1) specifying which column of data
; should appear in the plot.
;
  z_factor = 0.5
  z_column = fix(nz*z_factor)
  ez_line = ez(z_column,*)

  print, ' '
  print, 'z_column = ', z_column
  help, ez_line

;
; calculate the 10% level
;
  ez_lineMax = max(ez_line)
  ez_lineMin = ez_lineMax/10000.
  ez_level   = ez_lineMax/10.
  whez   = where( ez_line gt ez_level )
  numHighInd = n_elements(whez)
  y1 = ez_line(whez(numHighInd-1)+1)
  x1 = r(whez(numHighInd-1)+1)
  y2 = ez_line(whez(numHighInd-1))
  x2 =  r(whez(numHighInd-1))
  a0 = (x2*y1 - x1*y2)/(x2 - x1)
  a1 = (y2 - y1)/(x2 - x1)
  xlevel = (ez_level - a0)/a1
  print, 'Blue shift = ', xlevel

; Put the plot into a new window
  window_number = !d.window + 1
  print, ' '
  print, 'Line plot will appear in window ', window_number
  window, window_number
  lineout_i = 0
  lineout_jump:

; Load in the STD GAMMA-II color table
  loadct, 5

; Specify a font that looks great for printing (crappy on screen),
;   or else one that looks OK on the screen (also OK for printing).
  if (!d.name eq 'PS') then begin
    !p.font=1
    !p.charsize=1.6
    !p.charthick=1.5
    print, ' '
    print, 'Writing the 2D lineout plot to file ' + ps_line_file
  endif else begin
    !p.font=-1
    !p.charsize=2.0
    device,decomposed=0
    print, ' '
    print, 'Rendering the 2D lineout plot for the screen...'
  endelse

if (log_plot_flag eq 0) then begin
  ;plot, r, ez_line, thick=2, /xs, xtitle=y_label, ytitle=z_label, yrange=[ez_lineMin, ez_lineMax]
  plot, r, ez_line, thick=2, /xs, xtitle=y_label, ytitle=z_label, xrange=[waveLenMin, waveLenMax], yrange=[ez_lineMin, ez_lineMax]
  PlotS, [xlevel, xlevel], [ez_lineMin, ez_lineMax], lineStyle = 1
  PlotS, [xlevel, xlevel], [ez_lineMin, ez_lineMax], lineStyle = 1
endif else begin
  ;plot_io, r, ez_line, thick=2, /xs, xtitle=y_label, ytitle=z_label, yrange=[ez_lineMin, ez_lineMax]
  plot_io, r, ez_line, thick=2, /xs, xtitle=y_label, ytitle=z_label, xrange=[waveLenMin, waveLenMax], yrange=[ez_lineMin, ez_lineMax]
  PlotS, [waveLenMin, waveLenMax], [ez_level, ez_level], lineStyle = 1
  PlotS, [xlevel, xlevel], [ez_lineMin, ez_lineMax], lineStyle = 1 
  ;strout = '10 % level = ' + Fix(xlevel)
  ;XYOutS, 0.3, 0.7, strout, Size=2., /Normal
  XYOutS, 0.3,  0.7, '10 % level = ', Size=2., /Normal
  XYOutS, 0.46, 0.7, Fix(xlevel),  Size=2., /Normal
  XYOutS, 0.57, 0.7, '(nm)', Size=2., /Normal
endelse

; rendering to screen is done.  here we render to ps file
  if (lineout_i eq 0) then begin
    lineout_i = 1
    set_plot, 'PS', /copy
    device, filename = ps_line_file, /color, bits_per_pixel=8
    goto, lineout_jump  
  endif

; here rendering to file is done.  set back to screen defaults
  if (lineout_i eq 1) then begin
    device, /close
    set_plot, 'X'
    !p.font=-1
    !p.charsize=2.0
    device,decomposed=0
  endif

endif

; *****************************************************************
; Do a line plot of the inegrated along "y" power spectral density
; *****************************************************************

if (doIntegral eq 1) then begin
  nzFraction = 0.45  
  nzMin = fix(nzFraction*nz)
  nzMax = fix((1.-nzFraction)*nz) - 1
  if ( (nzMax + 1) eq nzMin ) then begin
    nzMax = nzMin
  endif 
  ;
  ; sum up along "y" 
  ; 
  psd_1d = dblarr(nr)
  for i = 0, nr-1 do begin
    psd_1d(i) = 0.
    for j = nzMin, nzMax do begin
      psd_1d(i) = psd_1d(i) + ez(j,i)
    endfor
  endfor

;  max_tmp = max(psd_1d)
;  psd_1d = psd_1d * max(ez) / max_tmp
;  ez_line = psd_1d 
  ez_line = psd_1d/(nzMax - nzMin + 1)

;
; calculate the 10% level
;
  ez_lineMax = max(ez_line)
  ez_lineMin = ez_lineMax/10000.
  ez_level   = ez_lineMax/10.
  whez   = where( ez_line gt ez_level )
  numHighInd = n_elements(whez)
  print, 'whez(numHighInd-1)+1 = ', whez(numHighInd-1)+1
  y1 = ez_line(whez(numHighInd-1)+1)
  x1 = r(whez(numHighInd-1)+1)
  y2 = ez_line(whez(numHighInd-1))
  x2 =  r(whez(numHighInd-1))
  a0 = (x2*y1 - x1*y2)/(x2 - x1)
  a1 = (y2 - y1)/(x2 - x1)
  xlevel = (ez_level - a0)/a1
  print, 'Blue shift = ', xlevel

; Put the plot into a new window
  window_number = !d.window + 1
  print, ' '
  print, 'Line plot will appear in window ', window_number
  window, window_number
  lineout_i = 0
  integral_jump:

; Load in the STD GAMMA-II color table
  loadct, 5

; Specify a font that looks great for printing (crappy on screen),
;   or else one that looks OK on the screen (also OK for printing).
  if (!d.name eq 'PS') then begin
    !p.font=1
    !p.charsize=1.6
    !p.charthick=1.5
    print, ' '
    print, 'Writing the 2D lineout plot to file ' + ps_line_file
  endif else begin
    !p.font=-1
    !p.charsize=2.0
    device,decomposed=0
    print, ' '
    print, 'Rendering the 2D lineout plot for the screen...'
  endelse

if (log_plot_flag eq 0) then begin
  ;plot, r, ez_line, thick=2, /xs, xtitle=y_label, ytitle=z_label, yrange=[ez_lineMin, ez_lineMax]
  plot, r, ez_line, thick=2, /xs, xtitle=y_label, ytitle=z_label, $
    xrange=[waveLenMin, waveLenMax], yrange=[ez_lineMin, ez_lineMax], $
    Title='Averaged 1d Power Spectral Density'
  PlotS, [xlevel, xlevel], [ez_lineMin, ez_lineMax], lineStyle = 1
  PlotS, [xlevel, xlevel], [ez_lineMin, ez_lineMax], lineStyle = 1
endif else begin
  ;plot_io, r, ez_line, thick=2, /xs, xtitle=y_label, ytitle=z_label, yrange=[ez_lineMin, ez_lineMax]
  plot_io, r, ez_line, thick=2, /xs, xtitle=y_label, ytitle=z_label, $
    xrange=[waveLenMin, waveLenMax], yrange=[ez_lineMin, ez_lineMax], $
    Title='Averaged 1d Power Spectral Density'
  PlotS, [waveLenMin, waveLenMax], [ez_level, ez_level], lineStyle = 1
  PlotS, [xlevel, xlevel], [ez_lineMin, ez_lineMax], lineStyle = 1 
  ;strout = '10 % level = ' + Fix(xlevel)
  ;XYOutS, 0.3, 0.7, strout, Size=2., /Normal
  XYOutS, 0.3,  0.7, '10 % level = ', Size=2., /Normal
  XYOutS, 0.46, 0.7, Fix(xlevel),  Size=2., /Normal
  XYOutS, 0.57, 0.7, '(nm)', Size=2., /Normal
endelse

; rendering to screen is done.  here we render to ps file
  if (lineout_i eq 0) then begin
    lineout_i = 1
    set_plot, 'PS', /copy
    device, filename = ps_integral_file, /color, bits_per_pixel=8
    goto, integral_jump  
  endif

; here rendering to file is done.  set back to screen defaults
  if (lineout_i eq 1) then begin
    device, /close
    set_plot, 'X'
    !p.font=-1
    !p.charsize=2.0
    device,decomposed=0
  endif

endif

; *****************************************************************
; Final clean-up
; *****************************************************************

; Stop here so that all variables will still be in scope.
; stop

; All done.
end





; *****************************************************************
; checkComponent procedure
; *****************************************************************

pro checkComponent, componentLabel, component

  CASE ComponentLabel OF
    ; if componentLabel is 'x' 
    'x': print, "Setting "+component+" label to "+ComponentLabel

    ; 'y'
    'y': print, "Setting "+component+" label to "+ComponentLabel

    ; 'r'
    'r': print, "Setting "+component+" label to "+ComponentLabel
                         
    ; 'z'                
    'z': print, "Setting "+component+" label to "+ComponentLabel

    ; no matches => print warning
  ELSE: print, "WARNING: Component label is not in the set (x, y, r, z)!?"
  ENDCASE
end