plot_postscript.f90

Spectral Element Method for wave propagation and rupture dynamics.

! SEM2DPACK version 2.2.11 -- A Spectral Element Method for 2D wave propagation and fracture dynamics,
!                             with emphasis on computational seismology and earthquake source dynamics.
! 
! Copyright (C) 2003-2007 Jean-Paul Ampuero
! All Rights Reserved
! 
! Jean-Paul Ampuero
! 
! ETH Zurich (Swiss Federal Institute of Technology)
! Institute of Geophysics
! Seismology and Geodynamics Group
! ETH H鰊ggerberg HPP O 13.1
! CH-8093 Z黵ich
! Switzerland
! 
! [e-mail]
! [phone] (office)
! [phone] (fax)
! 
! http://www.sg.geophys.ethz.ch/geodynamics/ampuero/
! 
! 
! This software is freely available for scientific research purposes. 
! If you use this software in writing scientific papers include proper 
! attributions to its author, Jean-Paul Ampuero.
! 
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License
! as published by the Free Software Foundation; either version 2
! of the License, or (at your option) any later version.
! 
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
! GNU General Public License for more details.
! 
! You should have received a copy of the GNU General Public License
! along with this program; if not, write to the Free Software
! Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
! 
module plot_postscript

  use spec_grid 
  use constants, only : NDIME

  implicit none
  private

  double precision, parameter :: centim = 28.5d0
  integer, parameter :: maxcolors = 20
  real, dimension(3,maxcolors), save :: RGB
  logical, save :: legend,interpol,vectors,mesh,symbols,boundaries,color,numbers
  integer, save :: set_background

  double precision, allocatable, save :: se_interp(:,:,:), fe_interp(:,:,:)

  integer, save :: isubsamp,DisplayPts
                   
  logical, parameter :: usletter=.true. ! Page format: US letter or A4

  double precision, save :: sizex,sizez,rapp_page,xmin,zmin,xmax,zmax,ScaleField
  character(4), save :: version='2.11'

  real, save :: usoffset

  public :: PLOT_PS,POST_PS_read,POST_PS_init

contains

!=======================================================================
!
!
! BEGIN INPUT BLOCK
!
! NAME   : PLOTS_POSTCRIPT
! GROUP  : PLOTS
! PURPOSE: Preferences for PostScript snapshots
! SYNTAX : &PLOTS_POSTSCRIPT vectors, mesh, background, color,
!               isubsamp, boundaries, symbols, numbers, legend,
!               ScaleField, Interpol, DisplayPts /
!
! ARG: vectors          [log] [F] Plots a vectorial field with arrows
! ARG: mesh             [log] [F] Plots the mesh on background
! ARG: background       [char] [''] Filled background, only for vector plots:
!                                   ''   none 
!                                   'P'  P-velocity model
!                                   'S'  S-velocity model
!                                   'T'  domains 
! ARG: isubsamp         [int] [3] Subsampling of the GLL nodes for the
!                                 output of velocity model. 
!                                 The default samples every 3 GLL points.
! ARG: boundaries       [log] [T] Colors every tagged boundary
! ARG: symbols          [log] [T] Plots symbols for sources and receivers
! ARG: numbers          [log] [F] Plots the element numbers
! ARG: legend           [log] [T] Writes legends
! ARG: color            [log] [T] Color output
! ARG: ScaleField       [dble] [0d0] Fixed amplitude scale (saturation),
!                       convenient for comparing snapshots and making movies. 
!                       The default scales each snapshot by its maximum amplitude
! ARG: Interpol         [log] [T] Interpolate field on a regular subgrid 
!                       inside each element
! ARG: DisplayPts       [log] [3] Size of interpolation subgrid inside each 
!                       element is DisplayPts*DisplayPts. The default plots at 
!                       vertices, mid-edges and element center.
!               
! END INPUT BLOCK

  subroutine POST_PS_read(iin)

  use echo, only : iout,echo_input

  integer, intent(in) :: iin
  character :: background
  character(10) :: bg_name

  NAMELIST / PLOTS_POSTSCRIPT / vectors,numbers &
                         ,background,isubsamp,color &
                         ,boundaries,symbols,legend,mesh &
                         ,interpol,DisplayPts,ScaleField

  vectors    = .false. 
  numbers    = .false. 
  background = ''
  isubsamp   = 3
  color      = .true.
  boundaries = .true.
  symbols    = .true.
  legend     = .true.
  mesh       = .false.
  interpol   = .true.
  DisplayPts = 3
  ScaleField = 0.d0 
     
  rewind(iin)
  read(iin,PLOTS_POSTSCRIPT,END=100)
100 continue
 
  if (vectors) then
    select case(background)
      case('P'); set_background=1; bg_name='P model'
      case('S'); set_background=2; bg_name='S model'
      case('T'); set_background=3; bg_name='domains'
      case default;   set_background=0; bg_name='none'
    end select
  else
    set_background=0
    bg_name = 'none'
  endif

  if (echo_input) write(iout,200) mesh,numbers &
                         ,bg_name,isubsamp,color &
                         ,boundaries,symbols,legend &
                         ,vectors,ScaleField,interpol,DisplayPts
  return

  200 format(//1x,'P o s t S c r i p t   O u t p u t s',/1x,35('='),//5x, &
  'Plot mesh . . . . . . . . . . . . . . . . . . (mesh) = ',L1/ 5x, &
  'Numbered mesh . . . . . . . . . . . . . . .(numbers) = ',L1/ 5x, &
  'Background fill . . . . . . . . . . . . (background) = ',A/ 5x, &
  'Subsampling for velocity model display  . (isubsamp) = ',I0/5x, &
  'Color display . . . . . . . . . . . . . . . .(color) = ',L1/ 5x, &
  'Plot boundaries . . . . . . . . . . . . (boundaries) = ',L1/ 5x, &
  'Plot symbols  . . . . . . . . . . . . . . .(symbols) = ',L1/ 5x, &
  'Write legends . . . . . . . . . . . . . . . (legend) = ',L1/ 5x, &
  'Plot vector fields  . . . . . . . . . . . .(vectors) = ',L1/ 5x, &
  'Amplitude-Scaling . . . . . . . . . . . (ScaleField) = ',F0.2/5x, &
  'Interpolate vector field  . . . . . . . . (interpol) = ',L1/5x, &
  'Points per edge for interpolation . . . (DisplayPts) = ',I0)


  end subroutine POST_PS_read

!=======================================================================
! NOTE: indexed color scales can also be defined in Level 2 PostScript by
!       [/Indexed /DeviceRGB 255 <... ... ...> ] setcolorspace 

  subroutine POST_PS_init(grid)

  use spec_grid, only : sem_grid_type,SE_init_interpol
  use fem_grid, only : FE_getshape,Fe_GetNodesPerElement

  type(sem_grid_type), intent(in) :: grid

  !-- figure size as % of page
  double precision, parameter :: rpercentx = 70.0d0, rpercentz = 77.0d0

  double precision :: xi , eta
  integer :: i,j,ngnod

  !-- color palette

  RGB(:,1)  = (/ 1.00, 0.00, 0.00 /) ! red
  RGB(:,2)  = (/ 0.00, 0.00, 1.00 /) ! blue
  RGB(:,3)  = (/ 0.93, 0.51, 0.93 /) ! violet
  RGB(:,4)  = (/ 0.73, 0.33, 0.83 /) ! medium orchid
  RGB(:,5)  = (/ 0.60, 0.20, 0.80 /) ! dark orchid
  RGB(:,6)  = (/ 0.54, 0.17, 0.89 /) ! blue violet
  RGB(:,7)  = (/ 0.42, 0.35, 0.80 /) ! slate blue
  RGB(:,8)  = (/ 1.00, 0.08, 0.58 /) ! deep pink
  RGB(:,9)  = (/ 0.12, 0.56, 1.00 /) ! dodger blue
  RGB(:,10) = (/ 0.00, 0.81, 0.82 /) ! dark turquoise
  RGB(:,11) = (/ 0.25, 0.88, 0.82 /) ! turquoise
  RGB(:,12) = (/ 0.20, 0.80, 0.20 /) ! lime green
  RGB(:,13) = (/ 0.00, 1.00, 0.50 /) ! spring green
  RGB(:,14) = (/ 0.50, 1.00, 0.00 /) ! chartreuse
  RGB(:,15) = (/ 0.68, 1.00, 0.18 /) ! green yellow
  RGB(:,16) = (/ 1.00, 1.00, 0.00 /) ! yellow
  RGB(:,17) = (/ 1.00, 0.98, 0.80 /) ! lemon chiffon
  RGB(:,18) = (/ 1.00, 0.84, 0.00 /) ! gold
  RGB(:,19) = (/ 1.00, 0.89, 0.71 /) ! mocassin
  RGB(:,20) = (/ 1.00, 0.85, 0.73 /) ! peach puff

  ! paper format A4 or US letter
  if(usletter) then
    usoffset = 1.75
    sizex = 27.94d0
    sizez = 21.59d0
  else
    usoffset = 0.
    sizex = 29.7d0
    sizez = 21.d0
  endif

  !-- grid extrema
  xmax=maxval(grid%coord(1,:))
  xmin=minval(grid%coord(1,:))
  zmax=maxval(grid%coord(2,:))
  zmin=minval(grid%coord(2,:))

  ! page size / physical domain size
  rapp_page = min(rpercentz*sizez/(zmax-zmin),rpercentx*sizex/(xmax-xmin))/100.d0

 ! weigths for interpolation (regular subgrids on each element)
  if (interpol) then
    ngnod = FE_GetNodesPerElement(grid%fem)
    allocate(se_interp(grid%ngll*grid%ngll,DisplayPts,DisplayPts))
    allocate(fe_interp(ngnod,DisplayPts,DisplayPts))
    do j=1,DisplayPts
      eta = 2.d0* dble(j-1)/dble(DisplayPts-1) - 1.d0
      do i=1,DisplayPts
        xi = 2.d0* dble(i-1)/dble(DisplayPts-1) - 1.d0
        call SE_init_interpol(xi,eta,se_interp(:,i,j),grid) 
        fe_interp(:,i,j) = FE_getshape(xi,eta,ngnod)
      enddo
    enddo
  endif

 
  end subroutine POST_PS_init


!=======================================================================
! PostScript plot manager
!
! Cases:
!   vfield is present and is vectorial: 
!       if requested ('vectors') plot as arrows and
!       fill the elements with efield, velocity model, domain tags or none
!       else plot the amplitude as color cells
!   vfield is present and is scalar: color cells
!   efield is present: if vfield is present, use as background,
!       else plot as color cells
!
  subroutine PLOT_PS(file,vfield,efield,grid,elast,stitle &
                    ,it_in,time_in,src,rec,comp)

  use stdio, only : IO_new_unit
  use echo, only : echo_run,iout,fmt1,fmtok
  use elastic, only : elast_type
  use sources, only : source_type,SO_inquire
  use receivers, only : rec_type,REC_inquire
  use fem_grid, only : FE_GetNodesPerElement, FE_GetElementCoord

  character(*)       , intent(in) :: file 
  type(sem_grid_type), intent(in) :: grid
  double precision   , intent(in), optional :: vfield(:,:)
  double precision   , intent(in), optional :: efield(grid%nelem)
  type(elast_type)   , intent(in) :: elast
  character(*)       , intent(in) :: stitle
  integer            , intent(in), optional :: it_in
  double precision   , intent(in), optional :: time_in
  type(source_type)  , optional, pointer :: src(:)
  type(rec_type)     , optional, pointer :: rec
  integer            , optional :: comp

 ! height of the domain tags/numbers, in centimeters
  double precision, parameter :: height = 0.25d0 

  double precision :: maxfield,time
  integer :: psunit,it,opt_comp
  logical :: nodal_field,elem_field,vector_field

  if (present(time_in) .and. present(it_in)) then
    time = time_in
    it = it_in
  else
    time = 0d0
    it = 0
  endif

  if (present(comp)) then
    opt_comp = comp
  else
    opt_comp = 0
  endif

  nodal_field = present(vfield)
  if (nodal_field) then
    maxfield = maxval(abs(vfield))
    vector_field = (size(vfield,2)==2)
  endif

  elem_field = present(efield)
  if (elem_field) maxfield = maxval(abs(efield))

  psunit = IO_new_unit()
  open(unit=psunit,file=trim(file),status='unknown')
!  if (echo_run) write(iout,'("Dump PostScript ",A," ...")',advance='no') trim(file)
  if (echo_run) write(iout,fmt1,advance='no') "Dump PostScript "//trim(file)

  call plot_header()

  if (legend) call plot_legend()

  ! Scale factor in the X and Z direction
  write(psunit,*) '%'
  write(psunit,*) '1.0 1.0 scale'
  write(psunit,*) '%'

  if (nodal_field) then

    if (vector_field) then
     !-- plot a nodal vector field with arrows
      if (vectors) then
        if (elem_field) then
          call plot_efield()
        else 
          call plot_model(set_background)
        endif
        if (mesh) call plot_mesh(color,numbers, set_background==3 )
        call plot_boundaries()
        call plot_vect()
     !-- plot the amplitude of a nodal vector field
      else
        if (opt_comp<1 .or. opt_comp>2) then
          call plot_scal( sqrt(vfield(:,1)*vfield(:,1)+vfield(:,2)*vfield(:,2)) )
        else
          call plot_scal( vfield(:,opt_comp) )
        endif
        if (mesh) call plot_mesh(color,numbers, .false.)