plot_efield.pro

pic 模拟程序！面向对象

; NAME:      plot_efield
;
; PURPOSE:   IDL procedure to read in E-field surface data from a .txt file
;            generated by dumping ASCII data from the XOOPIC gui.
;
; CATEGORY:
;
; CALLING SEQUENCE:
;   plot_efield, f_name, x_name, y_name, z_name, ScaleFlag, ScaleInput
;
; 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")
;  ScaleFlag:  a flag to specify the scaling for length, i.e. of q:
;                  0 -> no scaling
;                  1 -> length scale = c/omega_0 (the laser frequency)
;                  2 -> length scale = 1/Kp      (the plasma wave vector)
;                  3 -> length scale = lambda_p  (the plasma wavelength)
;  ScaleInput; for ScaleFlag = 0, input dummy value
;              for ScaleFlag = 1, input omega_0
;              for ScaleFlag = 2, input plasma density in [1/cm^3]
;              for ScaleFlag = 3, input plasma density in [1/cm^3]
;
; OPTIONAL INPUTS:     None.
;
; KEYWORD PARAMETERS:  None.
;
; OUTPUTS:             None.
;
; OPTIONAL OUTPUTS:    None.
;
; COMMON BLOCKS:       None.
;
; SIDE EFFECTS:        None.
;
; RESTRICTIONS: 
;
;    The XOOPIC .txt files generally have some lines of text.  You
;      need to delete any such text lines, leaving only the data.
;
;    The argument base_name is a string that specifies the input
;      data file (base_name.txt) and is also used to name output files.
;
;    This procedure is specifically for use with the Ex data, but it
;      could be easily modified for any other type of surface data:
;      a) Rename the _file and _label variables appropriately.
;      b) Look for the "convert" comment below and change the scaling
;         of the data to something appropriate for the desired units.
;
; PROCEDURE:
;
; EXAMPLE:
;   plot_efield, 'ey_file', 'x', 'y', 'y', 2, 1.93e18
;
; MODIFICATION HISTORY:
;   Nov  7, 2000:  original code (DLB)
;   May  1, 2001:  some generalization (DAD)
;   May 10, 2001:  renamed, with many changes (DLB)
;   Nov  2, 2001:  minor changes (DLB)
;   Sep 23, 2003:  minor changes (DLB)
;
;   Version: $Id: plot_efield.pro,v 1.12 2003/11/07 22:06:15 bruhwile Exp $
;
;   Copyright (c) 2000-2001 by Tech-X Corporation.  All rights reserved.


pro plot_efield, base_name, x1Label, x2Label, EcomponentLabel, ScaleFlag, ScaleInput

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

; Specify the text, restore and postscript file names
text_file    = base_name + '.txt'
restore_file = base_name + '.dat'
ps_ccon_file = base_name + '_colorcon.ps'
ps_cont_file = base_name + '_contour.ps'
ps_surf_file = base_name + '_surface.ps'
ps_line_file = base_name + '_lineout.ps'
ps_tier_file = base_name + '_tiered.ps'
lineout_file = base_name + '_lineout.dat'

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

doFFT      = 0
doShow3D   = 0
doSurface  = 0
doLineout  = 1
doContour  = 0
doColorCon = 1

; Specify non-default window sizes for contour, surface, show3d plots
; Dimension is pixels for the window command
x_win_ccon = 1000
y_win_ccon =  350
; Dimension is cm for the device command
x_device = 25.0
y_device = 12.0

; *****************************************************************
; Specify minimum and/or maximum values for x1 and x2:
; *****************************************************************
; Set auto_scale_min to 1 (default choice) for automatic scaling
;    of x1min and x2min.
; Set auto_scale_max to 1 (default choice) for automatic scaling
;    of x1max and x2max.

; If you chose to set auto_scale_min or auto_scale_max to zero, then
;    you must provide the desired limits.

auto_scale_min = 0
auto_scale_max = 0

; Set a factor (between 0 and 1) specifying which row of data
;    (r=constant or y=constant) that you want for the lineout plot.
r_factor = 0.5

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

; Check to see if the "restore" file has been created:
spawn, "ls " + 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.
;  z1d =data.field1(2,*)
;  r1d =data.field1(3,*)
;  ez1d=data.field1(4,*)

  z1d =data.field1(0,*)
  r1d =data.field1(1,*)
  ez1d=data.field1(2,*)

  print, ' '
  print, 'z1d, r1d, ez1d are the columns of the raw data:'
  help,z1d
  help,r1d
  help,ez1d
  print, 'z1d(0) z1d(1) z1d(n-1) = ', z1d(0),z1d(1),z1d(n_elements(z1d)-1)
  print, 'r1d(0) r1d(1) r1d(n-1) = ', r1d(0),r1d(1),r1d(n_elements(r1d)-1)

; Extract the unique values for r and z grid points
  ztemp=z1d(sort(z1d))
  z=ztemp(uniq(ztemp))
  rtemp=r1d(sort(r1d))
  r=rtemp(uniq(rtemp))

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

  print, ' '
  print, 'z and r contain only the unique values of the original arrays:'
  help,z
  help,r
  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)

; Create a 2-D array that holds the gridded surface data
  print, ' '
  print, 'Creating ez (2D array) from the raw data (1D array)....'
  nx=nz
  ny=nr
  ez = dblarr(nx,ny)
  for i = 0, nx-1 do begin
    ez(i,*) = ez1d[i*ny:(i+1)*ny-1]
  endfor
  print, '  ...done! '
  help,ez

; Now rotate everything by 180 degrees to compare with Brad Shadwick's results
  ez = reverse(ez,1)

; Save so IDL doesn't have to repeatedly parse the ASCII file
  save, z,r,ez,nz,nr, filename = restore_file

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


; *****************************************************************
  if ( (EcomponentLabel EQ 'x') OR (EcomponentLabel EQ 'X') ) then begin
;   Now reverse the sign of Ex to compare with Brad Shadwick's results
    ez = -ez

; *****************************************************************
;   Now smooth the fields a bit
;    ez = smooth(ez, 3, /EDGE_TRUNCATE)
  endif


; *****************************************************************
; Scale the data so it corresponds to the desired units
; *****************************************************************
;
if ( ScaleFlag eq 1 ) then begin
  ;
  ; scale the length using as a c/omega_0 as a length scale
  ; 
  omega_0 = ScaleInput
  lengthScale = 3.0e8/omega_0
  E0 = 0.511e6 / lengthScale
  z = z/lengthScale
  r = r/lengthScale
  ez=ez/E0
  q1_label = '!3' + x1Label
  q2_label = '!3' + x2Label

  z_label = '!3E!D' + EcomponentLabel + '!N/E!D0!N'
  if ( (EcomponentLabel EQ 'x') OR (EcomponentLabel EQ 'X') ) then begin
    z_label = '!3' + '-E' + '!D' + EcomponentLabel + '!N/E!D0!N'
  endif
endif else if (ScaleFlag eq 2) then begin
  ;
  ; scale the length using 1/Kp, Kp is the plasma wave vector
  ; 
  density = ScaleInput
  Kp = 2*acos(-1.0)*9000.0*sqrt(density)/3.0e8
  Ep = 0.511e6 * Kp   
  z = z*Kp
  r = r*Kp
  ez=ez/Ep
  q1_label = '!3' + x1Label
  q2_label = '!3' + x2Label

  z_label = '!3E!D' + EcomponentLabel + '!N/E!Dp!N'
  if ( (EcomponentLabel EQ 'x') OR (EcomponentLabel EQ 'X') ) then begin
    z_label = '!3' + '-E' + '!D' + EcomponentLabel + '!N/E!Dp!N'
  endif
endif else if (ScaleFlag eq 3) then begin
  ;
  ; scale the length using lambda_p the plasma wavelength
  ; 
  density = ScaleInput
  lambda_p = 3.e8 / 9000.0 / sqrt(density)
  Kp = 2*acos(-1.0) / lambda_p
  Ep = 0.511e6 * Kp   
  z = z / lambda_p
  r = r / lambda_p
  ez=ez/Ep
  q1_label = '!3' + x1Label
  q2_label = '!3' + x2Label

  z_label = '!3E!D' + EcomponentLabel + '!N/E!Dp!N'
  if ( (EcomponentLabel EQ 'x') OR (EcomponentLabel EQ 'X') ) then begin
    z_label = '!3' + '-E' + '!D' + EcomponentLabel + '!N/E!Dp!N'
  endif
endif else begin
  ;
  ; this is the default case of no scaling 
  ;
  x_label = '!3' + x1Label + ' [m]'
  y_label = '!3' + x2Label + ' [m]'

  z_label = '!3E!D' + EcomponentLabel
  if ( (EcomponentLabel EQ 'x') OR (EcomponentLabel EQ 'X') ) then begin
    z_label = '!3' + '-E' + '!D' + EcomponentLabel
  endif
endelse

print, ' '
print, 'After normalization of the z and r arrays:'
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)

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

x1min_data = z(0)
x2min_data = r(0)
x1max_data = z(nz-1)
x2max_data = r(nr-1)

;x1min = 0.0
;x1max = z(nz-1) - 1.8151
x1min = 2.0
x1max = z(nz-1)

if ( (EcomponentLabel EQ 'z') OR (EcomponentLabel EQ 'Z') ) then begin
  x1min =  2.
  x1max = 12.
endif

x2min = 0.5 * r(nr-1) - 5
x2max = 0.5 * r(nr-1) + 5

help,x1min_data
help,x1max_data
help,x2min_data
help,x2max_data

if (auto_scale_min ne 0) then begin
  x1min = x1min_data
  x2min = x2min_data
endif

if (auto_scale_max ne 0) then begin
  x1max = x1max_data
  x2max = x2max_data
endif

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

help,x1min
help,x1max
help,x2min
help,x2max

; *****************************************************************
; 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(z, r, ez, x1min, x2min)
  z  = struct_A.xnew
  r  = struct_A.ynew
  ez = struct_A.znew

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

  print, ' '
  print, 'After applying the specified x1min and x2min --'
  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

if ( (x1max lt x1max_data) or (x2max lt x2max_data) ) then begin
  struct_A = array_cut_max(z, r, ez, x1max, x2max)
  z  = struct_A.xnew
  r  = struct_A.ynew
  ez = 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

; *****************************************************************
; Create a 2-D color map appropriate for the surface
; *****************************************************************

; Create a 2D array that specifies the color for each grid point
; The color map ranges from 1 to 254 (rather than 0 to 255) to
;    avoid some problems, which appear to be IDL bugs
;
ez_color = ( ez - min(ez) )
ez_color = 251. * ez_color / max(ez_color)

print, ' '
print, 'Color table ez_color:'
help, ez_color
print, 'min(ez_color) = ', min(ez_color)
print, 'max(ez_color) = ', max(ez_color)

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

if (doContour eq 1) then begin

; Get a new window
  window_number = !d.window + 1
  print, ' '
  print, 'Contour plot will appear in window ', window_number
  window, window_number, xsize=x_win_ccon, ysize=ywin_ccon, RETAIN=2
  contour_i = 0