diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES
index e93239e57a..72f98c217d 100644
--- a/src/DEPENDENCIES
+++ b/src/DEPENDENCIES
@@ -68,6 +68,7 @@ global.o: geometry_header.o
global.o: material_header.o
global.o: mesh_header.o
global.o: particle_header.o
+global.o: plot_header.o
global.o: source_header.o
global.o: tally_header.o
global.o: timing.o
@@ -106,6 +107,7 @@ input_xml.o: geometry_header.o
input_xml.o: global.o
input_xml.o: mesh_header.o
input_xml.o: output.o
+input_xml.o: plot_header.o
input_xml.o: string.o
input_xml.o: tally_header.o
input_xml.o: xml-fortran/templates/cross_sections_t.o
@@ -155,6 +157,7 @@ output.o: geometry_header.o
output.o: global.o
output.o: mesh_header.o
output.o: particle_header.o
+output.o: plot_header.o
output.o: string.o
output.o: tally_header.o
@@ -182,8 +185,12 @@ plot.o: error.o
plot.o: geometry.o
plot.o: geometry_header.o
plot.o: global.o
+plot.o: output.o
plot.o: particle_header.o
+plot.o: plot_header.o
+plot.o: ppmlib.o
plot.o: source.o
+plot.o: string.o
search.o: constants.o
search.o: error.o
diff --git a/src/Makefile b/src/Makefile
index 465845ff32..14b74bd85a 100644
--- a/src/Makefile
+++ b/src/Makefile
@@ -28,7 +28,7 @@ USE_HDF5 = no
ifeq ($(COMPILER),gnu)
F90 = gfortran
- F90FLAGS := -cpp -fbacktrace
+ F90FLAGS := -cpp -fbacktrace -DNO_F2008
LDFLAGS =
# Debugging
diff --git a/src/OBJECTS b/src/OBJECTS
index f817f4cd5b..e78fcd0c62 100644
--- a/src/OBJECTS
+++ b/src/OBJECTS
@@ -28,6 +28,8 @@ output.o \
particle_header.o \
physics.o \
plot.o \
+plot_header.o \
+ppmlib.o \
random_lcg.o \
search.o \
source.o \
@@ -35,4 +37,4 @@ source_header.o \
string.o \
tally.o \
tally_header.o \
-timing.o
\ No newline at end of file
+timing.o
diff --git a/src/geometry_header.F90 b/src/geometry_header.F90
index fbcf0353a9..aa1fcca5f6 100644
--- a/src/geometry_header.F90
+++ b/src/geometry_header.F90
@@ -64,6 +64,7 @@ module geometry_header
& surfaces(:) ! List of surfaces bounding cell -- note that
! parentheses, union, etc operators will be listed
! here too
+ integer, allocatable :: rgb(:) ! plotting color
end type Cell
! array index of universe 0
diff --git a/src/global.F90 b/src/global.F90
index 2404172f54..bd4b7471c4 100644
--- a/src/global.F90
+++ b/src/global.F90
@@ -9,6 +9,7 @@ module global
use material_header, only: Material
use mesh_header, only: StructuredMesh
use particle_header, only: Particle
+ use plot_header, only: Plot
use source_header, only: ExtSource
use tally_header, only: TallyObject, TallyMap, TallyScore
use timing, only: Timer
@@ -38,6 +39,7 @@ module global
type(Lattice), allocatable, target :: lattices(:)
type(Surface), allocatable, target :: surfaces(:)
type(Material), allocatable, target :: materials(:)
+ type(Plot), allocatable, target :: plots(:)
! Size of main arrays
integer :: n_cells ! # of cells
@@ -45,6 +47,7 @@ module global
integer :: n_lattices ! # of lattices
integer :: n_surfaces ! # of surfaces
integer :: n_materials ! # of materials
+ integer :: n_plots ! # of plots
! These dictionaries provide a fast lookup mechanism -- the key is the
! user-specified identifier and the value is the index in the corresponding
@@ -187,10 +190,6 @@ module global
! PLOTTING VARIABLES
logical :: plotting = .false.
- real(8) :: plot_origin(3)
- real(8) :: plot_width(2)
- real(8) :: plot_basis(6)
- real(8) :: pixel
! ============================================================================
! HDF5 VARIABLES
@@ -235,6 +234,7 @@ contains
if (allocated(lattices)) deallocate(lattices)
if (allocated(surfaces)) deallocate(surfaces)
if (allocated(materials)) deallocate(materials)
+ if (allocated(plots)) deallocate(plots)
! Deallocate cross section data, listings, and cache
if (allocated(nuclides)) deallocate(nuclides)
diff --git a/src/input_xml.F90 b/src/input_xml.F90
index 9d83c4b1df..e1a9222854 100644
--- a/src/input_xml.F90
+++ b/src/input_xml.F90
@@ -8,6 +8,8 @@ module input_xml
use global
use mesh_header, only: StructuredMesh
use output, only: write_message
+ use plot_header
+ use random_lcg, only: prn
use string, only: lower_case, to_str, str_to_int, str_to_real, &
split_string, starts_with, ends_with
use tally_header, only: TallyObject
@@ -259,6 +261,14 @@ contains
c % material = cell_(i) % material
c % fill = cell_(i) % fill
+ ! Set plot color
+ if (plotting) then
+ allocate(c % rgb(3))
+ c % rgb(1) = prn()*255
+ c % rgb(2) = prn()*255
+ c % rgb(3) = prn()*255
+ end if
+
! Check to make sure that either material or fill was specified
if (c % material == 0 .and. c % fill == 0) then
message = "Neither material nor fill was specified for cell " // &
@@ -529,6 +539,14 @@ contains
! Copy material id
m % id = material_(i) % id
+ ! Set plot color
+ if (plotting) then
+ allocate(m % rgb(3))
+ m % rgb(1) = prn()*255
+ m % rgb(2) = prn()*255
+ m % rgb(3) = prn()*255
+ end if
+
! Copy density -- the default value for the units is given in the
! material_t.xml file and doesn't need to be specified here, hence case
! default results in an error.
@@ -1086,16 +1104,17 @@ contains
end subroutine read_tallies_xml
!===============================================================================
-! READ_TALLIES_XML reads data from a tallies.xml file and parses it, checking
-! for errors and placing properly-formatted data in the right data structures
+! READ_PLOT_XML reads data from a plot.xml file
!===============================================================================
subroutine read_plot_xml
use xml_data_plot_t
- logical :: file_exists ! does tallies.xml file exist?
- character(MAX_LINE_LEN) :: filename ! absolute path to tallies.xml
+ integer i
+ logical :: file_exists ! does plot.xml file exist?
+ character(MAX_LINE_LEN) :: filename ! absolute path to plot.xml
+ type(Plot), pointer :: pl => null()
! Check if plot.xml exists
filename = trim(path_input) // "plot.xml"
@@ -1112,28 +1131,79 @@ contains
! Parse plot.xml file
call read_xml_file_plot_t(filename)
- ! Copy plotting origin
- if (size(origin_) == 3) then
- plot_origin = origin_
- end if
+ ! Allocate plots array
+ n_plots = size(plot_)
+ allocate(plots(n_plots))
- ! Copy plotting width
- if (size(width_) == 2) then
- plot_width = width_
- end if
+ do i = 1, n_plots
+ pl => plots(i)
- ! Read basis
- select case (basis_)
- case ("xy")
- plot_basis = (/ 1, 0, 0, 0, 1, 0 /)
- case ("yz")
- plot_basis = (/ 0, 1, 0, 0, 0, 1 /)
- case ("xz")
- plot_basis = (/ 1, 0, 0, 0, 0, 1 /)
- end select
+ ! Copy data into plots
+ pl % id = plot_(i) % id
+ pl % aspect = plot_(i) % aspect
- ! Read pixel width
- pixel = pixel_
+ if (size(plot_(i) % pixels) == 2) then
+ pl % pixels = plot_(i) % pixels
+ end if
+
+ select case (plot_(i) % color)
+ case ("cell")
+ pl % color = PLOT_COLOR_CELLS
+ case ("mat")
+ pl % color = PLOT_COLOR_MATS
+ case default
+ message = "Unsupported plot color '" // plot_(i) % color &
+ // "' in plot " // trim(to_str(i))
+ call fatal_error()
+ end select
+
+ select case (plot_(i) % type)
+ case ("slice")
+ pl % type = PLOT_TYPE_SLICE
+ !case ("points")
+ ! pl % type = PLOT_TYPE_POINTS
+ case default
+ message = "Unsupported plot type '" // plot_(i) % type &
+ // "' in plot " // trim(to_str(i))
+ call fatal_error()
+ end select
+
+ select case (plot_(i) % basis)
+ case ("xy")
+ pl % basis = PLOT_BASIS_XY
+ case ("xz")
+ pl % basis = PLOT_BASIS_XZ
+ case ("yz")
+ pl % basis = PLOT_BASIS_YZ
+ case default
+ message = "Unsupported plot basis '" // plot_(i) % basis &
+ // "' in plot " // trim(to_str(i))
+ call fatal_error()
+ end select
+
+ ! Copy plotting origin
+ if (size(plot_(i) % origin) == 3) then
+ pl % origin = plot_(i) % origin
+ else
+ message = "Origin must be length 3 " &
+ // "in plot " // trim(to_str(i))
+ call fatal_error()
+ end if
+
+ ! Copy plotting width
+ if (size(plot_(i) % width) == 3) then
+ pl % width = plot_(i) % width
+ else if (size(plot_(i) % width) == 2) then
+ pl % width(1) = plot_(i) % width(1)
+ pl % width(2) = plot_(i) % width(2)
+ else
+ message = "Bad plot width " &
+ // "in plot " // trim(to_str(i))
+ call fatal_error()
+ end if
+
+
+ end do
end subroutine read_plot_xml
diff --git a/src/main.F90 b/src/main.F90
index 4aa5633b9c..af9824a6a0 100644
--- a/src/main.F90
+++ b/src/main.F90
@@ -6,7 +6,7 @@ program main
use initialize, only: initialize_run
use intercycle, only: shannon_entropy, calculate_keff, synchronize_bank
use output, only: write_message, header
- use plot, only: run_plot
+ use plotter, only: run_plot
use physics, only: transport
use random_lcg, only: set_particle_seed
use source, only: get_source_particle
diff --git a/src/material_header.F90 b/src/material_header.F90
index 916578923b..d8e9ca4ea2 100644
--- a/src/material_header.F90
+++ b/src/material_header.F90
@@ -21,6 +21,8 @@ module material_header
character(12) :: sab_name ! name of S(a,b) table
integer :: sab_table = 0 ! index in sab_tables
integer :: sab_nuclide = 0 ! index of nuclide which has S(a,b) table
+
+ integer, allocatable :: rgb(:) ! plotting color
end type Material
end module material_header
diff --git a/src/output.F90 b/src/output.F90
index 0446ba5676..c4209ace13 100644
--- a/src/output.F90
+++ b/src/output.F90
@@ -10,6 +10,7 @@ module output
use global
use mesh_header, only: StructuredMesh
use particle_header, only: LocalCoord
+ use plot_header
use string, only: upper_case, to_str
use tally_header, only: TallyObject
@@ -891,21 +892,46 @@ contains
subroutine print_plot()
+ integer i
+ type(Plot), pointer :: pl => null()
+
! Display header for plotting
call header("PLOTTING SUMMARY")
- ! Print plotting origin
- write(ou,100) "Plotting Origin:", trim(to_str(plot_origin(1))) // &
- " " // trim(to_str(plot_origin(2))) // " " // &
- trim(to_str(plot_origin(3)))
+ do i=1,n_plots
+ pl => plots(i)
- ! Print plotting width
- write(ou,100) "Plotting Width:", trim(to_str(plot_width(1))) // &
- " " // trim(to_str(plot_width(2)))
+ ! Print plot id
+ write(ou,100) "Plot ID:", trim(to_str(pl % id))
- ! Print pixel width
- write(ou,100) "Pixel Width:", trim(to_str(pixel))
- write(ou,*)
+ ! Print plotting origin
+ write(ou,100) "Origin:", trim(to_str(pl % origin(1))) // &
+ " " // trim(to_str(pl % origin(2))) // " " // &
+ trim(to_str(pl % origin(3)))
+
+ ! Print plotting width
+ if (pl % type == PLOT_TYPE_SLICE) then
+
+ write(ou,100) "Width:", trim(to_str(pl % width(1))) // &
+ " " // trim(to_str(pl % width(2)))
+ write(ou,100) "Coloring:", trim(to_str(pl % color))
+ write(ou,100) "Basis:", trim(to_str(pl % basis))
+ write(ou,100) "Pixels:", trim(to_str(pl % pixels(1))) // " " // &
+ trim(to_str(pl % pixels(2)))
+
+ else if (pl % type == PLOT_TYPE_POINTS) then
+
+ write(ou,100) "Width:", trim(to_str(pl % width(1))) // &
+ " " // trim(to_str(pl % width(2))) // " " &
+ // trim(to_str(pl % width(3)))
+ write(ou,100) "Coloring:", trim(to_str(pl % color))
+ write(ou,100) "Ray Spacing:", trim(to_str(pl % aspect))
+
+ end if
+
+ write(ou,*)
+
+ end do
! Format descriptor for columns
100 format (1X,A,T25,A)
diff --git a/src/plot.F90 b/src/plot.F90
index ce7b26ed6c..e23ece1c1c 100644
--- a/src/plot.F90
+++ b/src/plot.F90
@@ -1,4 +1,4 @@
-module plot
+module plotter
use constants
use error, only: fatal_error
@@ -6,185 +6,172 @@ module plot
cross_lattice, cell_contains
use geometry_header, only: Universe, BASE_UNIVERSE
use global
+ use output, only: write_message
use particle_header, only: LocalCoord, deallocate_coord
+ use plot_header
+ use ppmlib, only: Image, init_image, allocate_image, &
+ deallocate_image, set_pixel
use source, only: initialize_particle
+ use string, only: to_str
implicit none
contains
+
!===============================================================================
-! RUN_PLOT generates a binary stream file containing a list of surface/lattice
-! crossings and what cell was traveled through. A Python script can then be used
-! to generate a plot based on the recorded crossings and cells
+! RUN_PLOT
!===============================================================================
subroutine run_plot()
- integer :: i ! loop index
- integer :: surface_crossed ! surface which particle is on
- integer :: lattice_crossed ! is surface crossing in lattice?
- integer :: last_cell ! most recent cell particle was in
- integer :: enter_surface ! entrance surface
- real(8) :: xyz(3) ! starting coordinates
- real(8) :: last_x_coord ! bounding x coordinate
- real(8) :: last_y_coord ! bounding y coordinate
- real(8) :: d ! distance to boundary
- real(8) :: distance ! distance particle travels
- logical :: found_cell ! found cell which particle is in?
- character(MAX_LINE_LEN) :: path_plot ! unit for binary plot file
- type(Cell), pointer :: c => null()
- type(Universe), pointer :: univ => null()
- type(LocalCoord), pointer :: coord => null()
+ integer :: i
+ type(Plot), pointer :: pl => null()
- ! Open plot file for binary writing
- path_plot = trim(path_input) // "plot.out"
- open(UNIT=UNIT_PLOT, FILE=path_plot, STATUS="replace", ACCESS="stream")
+ do i=1,n_plots
+ pl => plots(i)
- ! Write origin, width, basis, and pixel width to file
- write(UNIT=UNIT_PLOT) plot_origin
- write(UNIT=UNIT_PLOT) plot_width
- write(UNIT=UNIT_PLOT) plot_basis
- write(UNIT=UNIT_PLOT) pixel
+ ! Display output message
+ message = "Processing plot " // trim(to_str(pl % id)) // "..."
+ call write_message(5)
- ! Determine coordinates of the upper-left corner of the plot
- xyz(1) = plot_origin(1) - plot_width(1) / 2.0
- xyz(2) = plot_origin(2) + (plot_width(2) - pixel) / 2.0
- xyz(3) = plot_origin(3)
+ if (pl % type == PLOT_TYPE_SLICE) then
- ! Determine bounding x and y coordinates for plot
- last_x_coord = plot_origin(1) + plot_width(1) / 2.0
- last_y_coord = plot_origin(2) - plot_width(2) / 2.0
+ ! create 2d image
+ call create_ppm(pl)
+
+ end if
+ end do
+
+ end subroutine run_plot
+
+!===============================================================================
+! create_ppm
+!===============================================================================
+ subroutine create_ppm(pl)
+
+ type(Plot), pointer :: pl
+
+ type(Image) :: img
+ integer :: in_i, out_i
+ integer :: x, y
+ integer :: r, g, b
+ real(8) :: in_pixel, out_pixel
+ real(8) :: xyz(3)
+ logical :: found_cell
+ type(Cell), pointer :: c => null()
+
+ call init_image(img)
+ call allocate_image(img, pl % pixels(1), pl % pixels(2))
+
+ if (pl % basis == PLOT_BASIS_XY) then
+ in_i = 1
+ out_i = 2
+ in_pixel = pl % width(1)/dble(pl % pixels(1))
+ out_pixel = pl % width(2)/dble(pl % pixels(2))
+ xyz(1) = pl % origin(1) - pl % width(1) / 2.0
+ xyz(2) = pl % origin(2) - pl % width(2) / 2.0
+ xyz(3) = pl % origin(3)
+ else if (pl % basis == PLOT_BASIS_XZ) then
+ in_i = 1
+ out_i = 3
+ in_pixel = pl % width(1)/dble(pl % pixels(1))
+ out_pixel = pl % width(3)/dble(pl % pixels(2))
+ xyz(1) = pl % origin(1) - pl % width(1) / 2.0
+ xyz(2) = pl % origin(2)
+ xyz(3) = pl % origin(3) - pl % width(3) / 2.0
+ else if (pl % basis == PLOT_BASIS_YZ) then
+ in_i = 2
+ out_i = 3
+ in_pixel = pl % width(2)/dble(pl % pixels(1))
+ out_pixel = pl % width(3)/dble(pl % pixels(2))
+ xyz(1) = pl % origin(1)
+ xyz(2) = pl % origin(2) - pl % width(2) / 2.0
+ xyz(3) = pl % origin(3) - pl % width(3) / 2.0
+ end if
! allocate and initialize particle
allocate(p)
+ call initialize_particle()
+ p % coord % xyz = xyz
+ p % coord % uvw = (/ 1, 0, 0 /)
+ p % coord % universe = BASE_UNIVERSE
- ! loop over horizontal rays
- do while(xyz(2) > last_y_coord)
+ do y=1, img%height
+ do x=1, img%width
- ! initialize the particle and set starting coordinate and direction
- call initialize_particle()
+ call deallocate_coord(p % coord0 % next)
+ p % coord => p % coord0
- p % coord % xyz = xyz
- p % coord % uvw = (/ 1, 0, 0 /)
+ call find_cell(found_cell)
- ! write starting coordinate to file
- write(UNIT=UNIT_PLOT) p % coord % xyz
-
- ! Find cell that particle is currently in
- call find_cell(found_cell)
-
- ! =======================================================================
- ! MOVE PARTICLE FORWARD TO NEXT CELL
-
- if (.not. found_cell) then
- ! Clear any coordinates beyond first level
- call deallocate_coord(p % coord0 % next)
- p % coord => p % coord0
-
- distance = INFINITY
- univ => universes(BASE_UNIVERSE)
- do i = 1, univ % n_cells
- p % coord0 % xyz = xyz
- p % coord0 % cell = univ % cells(i)
-
- call distance_to_boundary(d, surface_crossed, lattice_crossed)
- if (d < distance) then
- ! Check to make sure particle is actually going into this cell
- ! by moving it slightly forward and seeing if the cell contains
- ! that coordinate
-
- p % coord0 % xyz = p % coord0 % xyz + (d + TINY_BIT) * p % coord0 % uvw
-
- c => cells(p % coord0 % cell)
- if (.not. cell_contains(c)) cycle
-
- ! Set new distance and retain pointer to this cell
- distance = d
- enter_surface = surface_crossed
- end if
- end do
-
- ! No cell was found on this horizontal ray
- if (distance == INFINITY) then
- p % coord0 % xyz(1) = last_x_coord
- write(UNIT_PLOT) p % coord0 % xyz, 0
-
- ! Move to next horizontal ray
- xyz(2) = xyz(2) - pixel
- cycle
- end if
-
- ! Write coordinate where next cell begins
- p % coord0 % xyz = xyz + distance * p % coord0 % uvw
- write(UNIT=UNIT_PLOT) p % coord0 % xyz, 0
-
- ! Process surface crossing for next cell
- p % coord0 % cell = NONE
- p % surface = -enter_surface
- call cross_surface(enter_surface)
- end if
-
- ! =======================================================================
- ! MOVE PARTICLE ACROSS HORIZONTAL TRACK
-
- do while (p % alive)
- ! save particle's current cell
- last_cell = p % coord % cell
-
- ! Calculate distance to next boundary
- call distance_to_boundary(distance, surface_crossed, lattice_crossed)
-
- ! Advance particle
- coord => p % coord0
- do while (associated(coord))
- coord % xyz = coord % xyz + distance * coord % uvw
- coord => coord % next
- end do
-
- ! If next boundary crossing is out of range of the plot, only include
- ! the visible portion and move to next horizontal ray
- if (p % coord0 % xyz(1) >= last_x_coord) then
- p % alive = .false.
- p % coord0 % xyz(1) = last_x_coord
-
- ! If there is no cell beyond this boundary, mark it as cell 0
- if (distance == INFINITY) p % coord % cell = 0
-
- ! Write ending coordinates to file
- write(UNIT=UNIT_PLOT) p % coord0 % xyz, last_cell
- cycle
- end if
-
- ! Write boundary crossing coordinates to file
- write(UNIT=UNIT_PLOT) p % coord0 % xyz, last_cell
-
- p % coord % cell = 0
- if (lattice_crossed /= NONE) then
- p % surface = NONE
- call cross_lattice(lattice_crossed)
+ if (.not. found_cell) then
+ r = 0
+ g = 0
+ b = 0
+ else
+ c => cells(p % coord % cell)
+ if (pl % color == PLOT_COLOR_MATS) then
+ r = materials(c % material) % rgb(1)
+ g = materials(c % material) % rgb(2)
+ b = materials(c % material) % rgb(3)
+ else if (pl % color == PLOT_COLOR_CELLS) then
+ r = c % rgb(1)
+ g = c % rgb(2)
+ b = c % rgb(3)
else
- p % surface = surface_crossed
- call cross_surface(last_cell)
-
- ! Since boundary conditions are disabled in plotting mode, we need
- ! to manually add the last segment
- if (surfaces(abs(surface_crossed)) % bc /= BC_TRANSMIT) then
- p % coord0 % xyz(1) = last_x_coord
- write(UNIT=UNIT_PLOT) p % coord0 % xyz, 0
- exit
- end if
+ r = 0
+ g = 0
+ b = 0
end if
+ end if
+ call set_pixel(img, x, y, r, g, b)
+
+ p % coord0 % xyz(in_i) = p % coord0 % xyz(in_i) + in_pixel
+ end do
+
+ p % coord0 % xyz(in_i) = xyz(in_i)
+ p % coord0 % xyz(out_i) = p % coord0 % xyz(out_i) + out_pixel
+ end do
+
+ call output_ppm(pl,img)
+
+ call deallocate_image(img)
+
+ end subroutine create_ppm
+
+
+!===============================================================================
+! output_ppm
+!===============================================================================
+ subroutine output_ppm(pl,img)
+
+ type(Plot), pointer :: pl
+ type(Image), intent(in) :: img
+
+ integer :: i, j
+ character(MAX_LINE_LEN) :: path_plot ! unit for binary plot file
+
+ path_plot = trim(path_input) // "slice" // trim(to_str(pl % id)) // ".ppm"
+ open(UNIT=UNIT_PLOT, FILE=path_plot)
+
+ write(UNIT_PLOT, '(A2)') 'P6'
+ write(UNIT_PLOT, '(I0,'' '',I0)') img%width, img%height
+ write(UNIT_PLOT, '(A)') '255'
+
+ do j=1, img%height
+ do i=1, img%width
+ write(UNIT_PLOT, '(3A1)', advance='no') achar(img%red(i,j)), &
+ achar(img%green(i,j)), &
+ achar(img%blue(i,j))
end do
-
- ! Move y-coordinate to next position
- xyz(2) = xyz(2) - pixel
end do
! Close plot file
close(UNIT=UNIT_PLOT)
+
+ end subroutine output_ppm
- end subroutine run_plot
-end module plot
+end module plotter
diff --git a/src/plot_header.F90 b/src/plot_header.F90
new file mode 100644
index 0000000000..bfe6f9b071
--- /dev/null
+++ b/src/plot_header.F90
@@ -0,0 +1,31 @@
+module plot_header
+
+ implicit none
+
+!===============================================================================
+! PLOT hold plot information
+!===============================================================================
+
+ type Plot
+ integer :: id ! Unique ID
+ integer :: type ! Type
+ integer :: color ! quantity to color regions by
+ real(8) :: origin(3) ! xyz center of plot location
+ real(8) :: aspect ! spacing between rays in raytracer
+ real(8) :: width(3) ! xyz widths of plot
+ integer :: basis ! direction of plot slice
+ integer :: pixels(2) ! pixel width/height of plot slice
+ end type Plot
+
+ integer :: PLOT_TYPE_SLICE = 1
+ integer :: PLOT_TYPE_POINTS = 2
+
+ integer :: PLOT_BASIS_XY = 1
+ integer :: PLOT_BASIS_XZ = 2
+ integer :: PLOT_BASIS_YZ = 3
+
+ integer :: PLOT_COLOR_CELLS = 1
+ integer :: PLOT_COLOR_MATS = 2
+
+
+end module plot_header
diff --git a/src/ppmlib.F90 b/src/ppmlib.F90
new file mode 100644
index 0000000000..2c9687617a
--- /dev/null
+++ b/src/ppmlib.F90
@@ -0,0 +1,101 @@
+module ppmlib
+
+ implicit none
+
+!===============================================================================
+! Image holds RGB information for output PPM image
+!===============================================================================
+
+ type Image
+ integer, dimension(:,:), pointer :: red, green, blue
+ integer :: width, height
+ end type Image
+
+contains
+
+ subroutine init_image(img)
+
+ type(Image), intent(out) :: img
+
+ nullify(img % red)
+ nullify(img % green)
+ nullify(img % blue)
+
+ img % width = 0
+ img % height = 0
+
+ end subroutine init_image
+
+ subroutine allocate_image(img, w, h)
+
+ type(Image), intent(inout) :: img
+ integer, intent(in) :: w, h
+
+ allocate(img % red(w, h))
+ allocate(img % green(w, h))
+ allocate(img % blue(w, h))
+
+ img % width = w
+ img % height = h
+
+ end subroutine allocate_image
+
+ subroutine deallocate_image(img)
+
+ type(Image) :: img
+
+ if ( associated(img % red) ) deallocate(img % red)
+ if ( associated(img % green) ) deallocate(img % green)
+ if ( associated(img % blue) ) deallocate(img % blue)
+ end subroutine deallocate_image
+
+
+ function inside_image(img, x, y) result(inside)
+
+ type(Image), intent(in) :: img
+ integer, intent(in) :: x, y
+ logical :: inside
+
+ inside = .false.
+
+ if ( ( x < img % width) .and. &
+ ( y < img % height ) .and. &
+ ( x >= 0 ) .and. &
+ ( y >= 0 ) ) inside = .true.
+
+ end function inside_image
+
+
+ function valid_image(img) result(valid)
+
+ type(Image), intent(in) :: img
+ logical :: valid
+
+ valid = .false.
+
+ if ( img % width == 0 ) return
+ if ( img % height == 0 ) return
+ if ( .not. associated(img % red) .or. &
+ .not. associated(img % green) .or. &
+ .not. associated(img % blue) ) return
+
+ valid = .true.
+
+ end function valid_image
+
+
+ subroutine set_pixel(img, x, y, r, g, b)
+
+ type(Image), intent(inout) :: img
+ integer, intent(in) :: x, y
+ integer, intent(in) :: r, g, b
+
+ if ( inside_image(img, x, y) .and. valid_image(img)) then
+ img % red(x+1,y+1) = mod(abs(r), 256)
+ img % green(x+1, y+1) = mod(abs(g), 256)
+ img % blue(x+1, y+1) = mod(abs(b), 256)
+ end if
+ end subroutine set_pixel
+
+
+end module ppmlib
diff --git a/src/xml-fortran/templates/plot_t.xml b/src/xml-fortran/templates/plot_t.xml
index b4dc9eaea8..ef0875c47a 100644
--- a/src/xml-fortran/templates/plot_t.xml
+++ b/src/xml-fortran/templates/plot_t.xml
@@ -3,9 +3,17 @@
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+