plot_fft_1d.pro

pic 模拟程序！面向对象

; NAME:      plot_FFT_1d
;
; PURPOSE:   IDL procedure to read in FFT power spectra density data from
;            a .txt file generated by dumping ASCII data from either the
;            XOOPIC or the OOPIC Pro GUI.
;
; CATEGORY:
;
; CALLING SEQUENCE:
;   plot_fft_1d, f_name, x_name, y_name, z_name
;
; INPUTS:
;   f_name:    string to specify base name of data file
;   x_name:    string to specify horizontal axis (must be "x" or "z")
;   y_name:    string to specify vertical   axis (must be "y" or "r")
;   z_name:    string to specify field component 
;              (must be "x", "y", "z", "r", or "phi")
;
; OPTIONAL INPUTS:     None.
;
; KEYWORD PARAMETERS:  None.
;
; OUTPUTS:             None.
;
; OPTIONAL OUTPUTS:    None.
;
; COMMON BLOCKS:       None.
;
; SIDE EFFECTS:        None.
;
; RESTRICTIONS: 
;
;    The XOOPIC and OOPIC Pro .txt files generally have some lines of text.
;    You need to delete any such text lines, leaving only the data.
;
;    XOOPIC and OOPIC Pro can generate both 1-D and 2-D FFT plots.
;    This procedure is for use with the 1-D FFT plots only.
;
; PROCEDURE:
;
; EXAMPLE:
;   plot_fft_1d, "ey_fft_file", "x", "y", "y"
;
; MODIFICATION HISTORY:
;   Oct 20, 2001:  initially adapted from the plot_FFT pro file (DLB)
;
;   Version: $Id: plot_FFT_1d.pro,v 1.1 2001/10/23 18:06:10 bruhwile Exp $
;
;   Copyright (c) 2001 by Tech-X Corporation.  All rights reserved.


pro plot_FFT_1d, base_name, x1Label, x2Label, EcomponentLabel

; *****************************************************************
; Set flag for log plots
; *****************************************************************

; If log_plot_flag is nonzero, then this procedure will take the
; logarithm of the z-data

log_plot_flag = 1


; *****************************************************************
; Names for files and labels
; *****************************************************************

; Specify the text, restore and postscript file names
text_file    = base_name + '.txt'
restore_file = base_name + '.dat'

if (log_plot_flag eq 0) then begin
  ps_ccon_file = base_name + '_colorcon_lin.ps'
  ps_cont_file = base_name + '_contour_lin.ps'
  ps_surf_file = base_name + '_surface_lin.ps'
  ps_line_file = base_name + '_integral_lin.ps'
  ps_tier_file = base_name + '_tiered_lin.ps'
endif else begin
  ps_ccon_file = base_name + '_colorcon_log.ps'
  ps_cont_file = base_name + '_contour_log.ps'
  ps_surf_file = base_name + '_surface_log.ps'
  ps_line_file = base_name + '_integral_log.ps'
  ps_tier_file = base_name + '_tiered_log.ps'
endelse

; *****************************************************************
; Specify which plots you want (1) or don't want (0):
; *****************************************************************

doIntegral = 1
doColorCon = 1
doContour = 0
doSurface = 0
doShow3D = 0

; *****************************************************************
; Specify minimum and/or maximum values for x1 and x2:
; *****************************************************************

; Set auto_scale_lambda to 1 for automatic scaling of the horizontal
;    axis (wavelength).
; lambda is in units of nm
;
auto_scale_lambda = 0
x1min = 469
x1max = 833

; Set auto_scale_y to 1 for automatic scaling of the y
;    axis (transverse position).
; y is in units of microns
;
auto_scale_y = 1
; x2min = 0.
; x2max = 200.

; Set auto_scale_count to 1 for automatic scaling of the vertical
;    axis ("count" or power spectral density).
; count is in "arbitrary units" (must be scaled)
;
auto_scale_count = 1
; x3min = 1.
; x3max = 13780.


; *****************************************************************
; Parse the ascii data file or restore from the binary IDL file
; *****************************************************************

; Check to see if the "restore" file has been created:
spawn, "ls | grep " + restore_file, check_file, /sh
;print, 'The "check"   file is: ', check_file
;print, 'The "restore" file is: ', restore_file

; If the restore file exists, then use it.
if (check_file(0) eq restore_file) then begin

  print, ' '
  print,'Reading from the restore file: ' + restore_file + ' ...'
  restore, filename = restore_file

; Otherwise, parse the text file:
endif else begin

  print, ' '
  print,'Parsing the text file: ' + text_file + ' ...'
  data=read_ascii(text_file)

  print, ' '
  print, 'Here is the size and shape of the raw data:'
  help,data.field1

; *****************************************************************
; Initial manipulation of the data
; *****************************************************************

; Load the raw data into individual arrays.
;
; Use the following for ascii data dumped from XOOPIC
;  y1d =data.field1(2,*)
;  k1d =data.field1(3,*)
;  temp=data.field1(4,*)
;
; Use the following for ascii data dumped from OOPIC Pro
  y1d =data.field1(2,*)
  k1d =data.field1(3,*)
  temp=data.field1(4,*)

  print, ' '
  print, 'y1d, k1d, temp are the columns of the raw data:'
  help,y1d
  help,k1d
  help,temp
  print, 'y1d(0) y1d(1) y1d(n-1) = ', y1d(0),y1d(1),y1d(n_elements(y1d)-1)
  print, 'k1d(0) k1d(1) k1d(n-1) = ', k1d(0),k1d(1),k1d(n_elements(k1d)-1)

; Extract the unique values for r and z grid points
  ytemp=y1d(sort(y1d))
  y=ytemp(uniq(ytemp))
  ktemp=k1d(sort(k1d))
  kx=ktemp(uniq(ktemp))

  ny = n_elements(y)
  nx = n_elements(kx)

; Create a 2-D array that holds the gridded surface data
  print, ' '
  print, 'Rearranging the structure of the raw data (can take a while)....'
  print, "max_temp = ", max(temp)
  print, "min_temp = ", min(temp)
  psd = dblarr(nx,ny)
  for i = 0, nx-1 do begin
     psd(i,*) = temp[i*ny:(i+1)*ny-1]
  endfor
  print, '  ...done! '

; strip out the kx(0) data, because this is not relevant,
;   and it leads to lambda(0)=inf below.
  kx = kx(1:nx-1)
  psd = psd(1:nx-1,*)
  nx = nx-1

  print, ' '
  print, 'nx = ', nx
  print, 'ny = ', ny
  print, ' '
  print, 'kx(0) kx(1) kx(', nx-1, ') = ', kx(0),kx(1),kx(nx-1)
  print, 'y(0)  y(1)  y( ', ny-1, ') = ', y(0), y(1), y(ny-1)
  print, ' '
  help,psd
  print, "max_psd = ", max(psd)
  print, "min_psd = ", min(psd)

; Save so IDL doesn't have to repeatedly parse the ASCII file
  save, y,kx,psd,nx,ny, filename = restore_file

; *****************************************************************
; Here is the end of the if/then/else construct for parsing.
; *****************************************************************
endelse


; *****************************************************************
; Define the axes and title labels
; *****************************************************************
;
z_label = '!3PSD(E!Dz!N)'
x_label = 'wavelength (nm)'
y_label = 'PSD (arbitrary units)'


; *****************************************************************
; Convert from wavenumber k to lambda
; *****************************************************************
;
; f(lambda)*d_lambda = f(k)*d_k ==> f(lambda) = f(k)*(d_k/d_lambda)
;

two_pi = 16.*acos(0.)
lambda = two_pi / kx
dk_dlambda = kx * kx / two_pi

for i = 0, nx-1 do begin
  psd(i,*) = psd(i,*) * dk_dlambda(i)
endfor

; *****************************************************************
; Normalize psd    to arbitrary units
; Normalize lambda from m to nm
; Normalize y      from m to microns
; *****************************************************************
  psdFactor = 1.0
  psd    = psdFactor * psd
  lambda = 1.0e+09   * lambda
  y      = 1.0e+06   * y

print, ' '
print, 'After normalization of the y, lambda and "count" arrays:'
print, ' '
print, 'nx = ', nx
print, 'ny = ', ny
print, ' '
print, 'lambda(0) lambda(1) lambda(', nx-1, ') = ', lambda(0),lambda(1),lambda(nx-1)
print, '     y(0)      y(1)      y(', ny-1, ') = ',      y(0),     y(1),     y(ny-1)
print, ' '
print, 'min_psd max_psd  = ', min(psd), max(psd)

print, ' '
print, 'These are the min and max values of the original data:'

x1min_data = lambda(0)
x1max_data = lambda(nx-1)
x2min_data = y(0)
x2max_data = y(ny-1)

if (auto_scale_lambda ne 0) then begin
  x1min = x1min_data
  x1max = x1max_data
endif

if (auto_scale_y ne 0) then begin
  x2min = x2min_data
  x2max = x2max_data
endif

; Sum up the integrated (over y) spectrum
;
psd_1d = dblarr(nx)
for i = 0, nx-1 do begin
  psd_1d(i) = 0.
  for j = 0, ny-1 do begin
    psd_1d(i) = psd_1d(i) + psd(i,j)
  endfor
endfor
max_tmp = max(psd_1d)
psd_1d = psd_1d * max(psd) / max_tmp

x3min_data = min(psd_1d)
x3max_data = max(psd_1d)

print, ' '
help, psd_1d
print, 'min_psd_1d = ', x3min_data
print, 'max_psd_1d = ', x3max_data

if (auto_scale_count ne 0) then begin
  x3min = x3min_data
  x3max = x3max_data
endif

print, ' '
print, 'These are the specified min and max values:'

help,x1min
help,x1max
help,x2min
help,x2max
help,x3min
help,x3max

; *****************************************************************
; Surface plots don't support xrange/yrange, so truncate the data:
; *****************************************************************
if ( (x1min gt x1min_data) or (x2min gt x2min_data) ) then begin
  struct_A = array_cut(lambda, y, psd, x1min, x2min)
  lambda = struct_A.xnew
  y      = struct_A.ynew
  psd    = struct_A.znew

  nx = n_elements(lambda)
  ny = n_elements(y)

  print, ' '
  print, 'After applying the specified x1min and x2min --'
  print, 'nx = ', nx
  print, 'ny = ', ny
  print, ' '
  print, 'lambda(0) lambda(1) lambda(', nx-1, ') = ', lambda(0),lambda(1),lambda(nx-1)
  print, '     y(0)      y(1)      y(', ny-1, ') = ',      y(0),     y(1),     y(ny-1)
endif

if ( (x1max lt x1max_data) or (x2max lt x2max_data) ) then begin
  struct_A = array_cut_max(lambda, y, psd, x1max, x2max)
  z  = struct_A.xnew
  r  = struct_A.ynew
  psd = struct_A.znew

  nz = n_elements(z)
  nr = n_elements(r)

  print, ' '
  print, 'After applying the specified x1max and x2max --'
  print, 'nz = ', nz
  print, 'nr = ', nr
  print, 'z(0) z(1) z(', nz-1, ') = ', z(0),z(1),z(nz-1)
  print, 'r(0) r(1) r(', nr-1, ') = ', r(0),r(1),r(nr-1)
endif


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

if (doIntegral eq 1) then begin

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

; Load in the STD GAMMA-II color table