#! /usr/bin/perl5
# Removed -w to get rid of "used only once" warnings

# In emacs, please make this -*- perl -*- mode. Thanks.

# Enhanced mcdisplay script for --trace output from McStas simulation.
#
# Implements graphic display using either
# PGPLOT
# Scilab
# Matlab
#
# PW 20030320
#
#   This file is part of the McStas neutron ray-trace simulation package
#   Copyright (C) 1997-2004, All rights reserved
#   Risoe National Laborartory, Roskilde, Denmark
#   Institut Laue Langevin, Grenoble, France
#
#   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






# Config module needed for Win32 vs unix setup.
# PW 20030320
use Config;

BEGIN {
  # default configuration (for all high level perl scripts)
  if($ENV{"MCSTAS"}) {
    $MCSTAS::sys_dir = $ENV{"MCSTAS"};
  } else {
    if ($Config{'osname'} eq 'MSWin32') {
      $MCSTAS::sys_dir = "c:\\mcstas\\lib";
    } else {
      $MCSTAS::sys_dir = "/usr/local/lib/mcstas";
    }
  }
  $MCSTAS::perl_dir = "$MCSTAS::sys_dir/tools/perl";

  # custom configuration (this script)

  # If this is Win32, load OLE related modules -
  # we can not talk to Matlab using pipe on Win32 :(
  # PW 20030314
  if ($Config{'osname'} eq 'MSWin32') {
    require Win32::OLE;
    require Win32::OLE::Variant;
  }
}

use lib $MCSTAS::perl_dir;
require "mcstas_config.perl";

require "mcrunlib.pl";
# IPC can probably be used safely, exists on sysv type systems,
# linux, Win32. Will investigate further regarding portability.
# PW 20030404
use IPC::Open2;

$magnification = 1;
$zooming = 0;

my (%transformations, @components);


sub read_instrument {
    my ($in) = @_;
    my ($st, $comp);

    $st = 0;
    @components = ();
    while(<$in>) {
        if($st == 0 && /^INSTRUMENT:/) {
            # Start of instrument description.
            $st = 1;
            if ($MCSTAS::mcstas_config{'PLOTTER'} =~ /Matlab/i) {
              # Initialize matlab struct...
              write_process("addpath('$MCSTAS::sys_dir/tools/matlab');\n");
              write_process("mcdisplay('Init');\n");
              write_process("global INSTRUMENT;\n");
              write_process("INSTRUMENT.descr='$sim';\n");
              # Possibly, set firstcomp + lastcomp
              if ($first) { write_process("INSTRUMENT.firstcomp='$first';\n"); }
              if ($last)  { write_process("INSTRUMENT.lastcomp ='$last';\n");  }
            }
            if ($MCSTAS::mcstas_config{'PLOTTER'} =~ /Scilab/i) {
              # Initialize scilab struct...
              write_process("exec('$MCSTAS::sys_dir/tools/scilab/mcdisplay.sci',-1);\n");
              write_process("INSTRUMENT.descr='$sim';\n");
              # Possibly, set firstcomp + lastcomp
              if ($first) { write_process("INSTRUMENT.firstcomp='$first';\n"); }
              if ($last)  { write_process("INSTRUMENT.lastcomp ='$last';\n");  }
              if ($save) {
                if ($direct_output) { write_process("INSTRUMENT.save_format='$direct_output';\n"); }
                else  { write_process("INSTRUMENT.save_format='-scg';\n"); }
                write_process("INSTRUMENT.save=1;\n");
              } else { write_process("INSTRUMENT.save=0;\n"); }
            }
        } elsif($st == 1 && /^COMPONENT:\s*"([a-zA-Z0-9_]+)"\s*/) {
            $comp = $1;
            @components = (@components, $comp);
            if ($MCSTAS::mcstas_config{'PLOTTER'} =~ /Matlab/i) {
              # Initialize components in matlab struct:
              write_process("INSTRUMENT.name{length(INSTRUMENT.name)+1}='$comp';\n");
            }
            if ($MCSTAS::mcstas_config{'PLOTTER'} =~ /Scilab/i) {
              # Initialize components in scilab struct:
              write_process("setcomponent('$comp');\n");
            }
        } elsif($st == 1 && /^POS:(.*)$/) {
            my @T;
            @T = split ",", $1;
            $transformations{$comp} = \@T;
            $compcnt=scalar(@components);
            if ($MCSTAS::mcstas_config{'PLOTTER'} =~ /Matlab/i) {
              # Component position for matlab struct:
              write_process("INSTRUMENT.$comp.T=ReshapeTransform([@T]);\n");
              write_process("INSTRUMENT.$comp.j=$compcnt;\n");
              write_process("INSTRUMENT.$comp.K=cell(0);\n");
            }
            if ($MCSTAS::mcstas_config{'PLOTTER'} =~ /Scilab/i) {
              # Component position for scilab struct:
              write_process("setposition([@T]);\n");
            }
        } elsif($st == 1 && /^MCDISPLAY: start$/) {
            $st = 2;                # Start of component graphics representation
        } elsif($st == 2 && /^MCDISPLAY: component ([a-zA-Z0-9_]+)/) {
            $comp = $1;
            $compdraw{$comp} = {};
            $compdraw{$comp}{'elems'} = [];
            if ($MCSTAS::mcstas_config{'PLOTTER'} =~ /Scilab/i) {
              # Initialize component variable etc. in scilab
              write_process("trace_component('$comp');\n");
            }
        } elsif($st == 2 && /^MCDISPLAY: magnify\('([xyz]*)'\)$/) {
            my $mag = $1;
            $compdraw{$comp}{'magX'} = 1 if $mag =~ /x/i;
            $compdraw{$comp}{'magY'} = 1 if $mag =~ /y/i;
            $compdraw{$comp}{'magZ'} = 1 if $mag =~ /z/i;
        } elsif($st == 2 && /^MCDISPLAY: multiline\(([0-9]+),([^()\n]+)\)$/) {
            my $count = $1;
            my @coords = split ',', $2;
            push @{$compdraw{$comp}{'elems'}},
                {type => 'multiline',
                 count => $count,
                 coords => \@coords};
            if ($MCSTAS::mcstas_config{'PLOTTER'} =~ /Matlab/i) {
              # Line elements for Matlab struct
              write_process("coords=[@coords];\n");
              write_process("x=coords(1:3:length(coords));\n");
              write_process("y=coords(2:3:length(coords));\n");
              write_process("z=coords(3:3:length(coords));\n");
              write_process("coords=[x;y;z;1+0*z];\n");
              write_process("INSTRUMENT.$comp.K{size(INSTRUMENT.$comp.K,2)+1}=coords;\n");
            }
            if ($MCSTAS::mcstas_config{'PLOTTER'} =~ /Scilab/i) {
              # Line elements for Scilab struct
              write_process("multiline([$1 @coords]);\n");
            }
        } elsif($st == 2 &&
                /^MCDISPLAY: circle\('([xyzXYZ]{2})',([-+0-9.eE]+),([-+0-9.eE]+),([-+0-9.eE]+),([-+0-9.eE]+)\)$/) {
            my ($plane,$x,$y,$z,$r) = ($1,$2,$3,$4,$5);
            # Make a circle using a 25-order multiline.
            my @coords = ();
            my $i;
            for($i = 0; $i <= 24; $i++) {
                my $a = $r*cos(2*3.1415927/24*$i);
                my $b = $r*sin(2*3.1415927/24*$i);
                my ($x1,$y1,$z1) = ($x,$y,$z);
                if($plane =~ /xy|yx/i) {
                    $x1 += $a;
                    $y1 += $b;
                } elsif($plane =~ /xz|zx/i) {
                    $x1 += $a;
                    $z1 += $b;
                } elsif($plane =~ /yz|zy/i) {
                    $y1 += $a;
                    $z1 += $b;
                } else {
                    die "mcdisplay: Bad plane specifier in circle: '$plane'";
                }
                push @coords, $x1, $y1, $z1;
            }
            push @{$compdraw{$comp}{'elems'}},
                {type => 'multiline',
                 count => 25,
                 coords => \@coords};
            if ($MCSTAS::mcstas_config{'PLOTTER'} =~ /Matlab/i) {
              # Line elements for Matlab struct, circle representation
              write_process("coords=[@coords];\n");
              write_process("x=coords(1:3:length(coords));\n");
              write_process("y=coords(2:3:length(coords));\n");
              write_process("z=coords(3:3:length(coords));\n");
              write_process("coords=[x;y;z;1+0*z];\n");
              write_process("INSTRUMENT.$comp.K{size(INSTRUMENT.$comp.K,2)+1}=coords;\n");
            }
            if ($MCSTAS::mcstas_config{'PLOTTER'} =~ /Scilab/i) {
              # Line elements for Scilab struct, circle representation
              write_process("circle('$plane',$x,$y,$z,$r);\n");
            }
        } elsif($st == 2 && /^MCDISPLAY: end$/) {
            $st = 1;  # End of component graphics representation
            if ($MCSTAS::mcstas_config{'PLOTTER'} =~ /Matlab/i) {
              # Matlab 'End of instrument'
              write_process("mcdisplay('Load');\n");
              write_process("PlotInstrument('init');\n");
              # Check if we were called with --save option, output matlab figure if so...
              if ($save) {
                # Clone the graph to another window...
                write_process("ax=gca;\n");
                write_process("h=figure('numbertitle','off','name','$sim McStas Instrument')\n;");
                write_process("copyobj(ax,h);\n");
                write_process("saveas(h,'$sim.fig','fig');\n");
                write_process("delete(h);\n");
                write_process("delete(INSTRUMENT.fig);\n");
                write_process("exit;\n");
              } else {
                write_process("set(INSTRUMENT.fig, 'closerequestfcn','exit;');\n");
                write_process("wait(INSTRUMENT.fig);\n");
              }
            }
            if ($MCSTAS::mcstas_config{'PLOTTER'} =~ /Scilab/i) {
              # Scilab 'End of instrument'
              write_process("endtrace();\n");
            }
        } elsif($st == 1 && /^INSTRUMENT END:/) {
            $st = 100;
            last;
        } else {
            print;
        }
    }
    exit if($st != 100);        # Stop when EOF seen before instrument end.
    return $#components + 1;
}


sub make_instrument {
    my (@x, @y, @z, @ori, @dis, @comp);
    my ($i, $c, %instr);
    my ($xmin, $xmax, $ymin, $ymax, $zmin, $zmax);

    $i = 0;
    foreach $c (@components) {
        my (@T, @U);
        @T = @{$transformations{$c}};
        $x[$i] = $T[0];
        $xmin = $x[$i] if !$xmin || $x[$i] < $xmin;
        $xmax = $x[$i] if !$xmax || $x[$i] > $xmax;
        $y[$i] = $T[1];
        $ymin = $y[$i] if !$ymin || $y[$i] < $ymin;
        $ymax = $y[$i] if !$ymax || $y[$i] > $ymax;
        $z[$i] = $T[2];
        $zmin = $z[$i] if !$zmin || $z[$i] < $zmin;
        $zmax = $z[$i] if !$zmax || $z[$i] > $zmax;

        @U = ($T[3], $T[4], $T[5], $T[6], $T[7], $T[8], $T[9], $T[10], $T[11]);
        $ori[$i] = \@U;
        $comp[$i] = $c;
        # Now transform coordinates for component graphics representations.
        if($compdraw{$c}) {
            my $magX = $compdraw{$c}{'magX'};
            my $magY = $compdraw{$c}{'magY'};
            my $magZ = $compdraw{$c}{'magZ'};
            foreach $elem (@{$compdraw{$c}{'elems'}}) {
                if($elem->{'type'} eq 'multiline') {
                    my @coords = @{$elem->{'coords'}};
                    my @xs = ();
                    my @ys = ();
                    my @zs = ();
                    my ($xv,$yv,$zv);
                    while(@coords) {
                        $xv = shift(@coords);
                        $yv = shift(@coords);
                        $zv = shift(@coords);
                        $xv *= $magnification if $magX;
                        $yv *= $magnification if $magY;
                        $zv *= $magnification if $magZ;
                        push @xs, ($xv*$T[3] + $yv*$T[6] + $zv*$T[9]  + $T[0]);
                        push @ys, ($xv*$T[4] + $yv*$T[7] + $zv*$T[10] + $T[1]);
                        push @zs, ($xv*$T[5] + $yv*$T[8] + $zv*$T[11] + $T[2]);
                    }
                    $elem->{'X'} = \@xs;
                    $elem->{'Y'} = \@ys;
                    $elem->{'Z'} = \@zs;
                }
            }
        }
        $i++;
    }
    if($xmin == $xmax) {
        $xmin--;
        $xmax++;
    }
    if($ymin == $ymax) {
        $ymin--;
        $ymax++;
    }
    if($zmin == $zmax) {
        $zmin--;
        $zmax++;
    }
    $xmin -= ($xmax - $xmin) / 6;
    $xmax += ($xmax - $xmin) / 6;
    $ymin -= ($xmax - $xmin) / 6;
    $ymax += ($xmax - $xmin) / 6;
    $zmin -= ($xmax - $xmin) / 6;
    $zmax += ($xmax - $xmin) / 6;
    %instr = ('x' => \@x, 'y' => \@y, z => \@z,
              ori => \@ori, dis => \@dis, comp => \@comp,
              xmin => $xmin, xmax => $xmax,
              ymin => $ymin, ymax => $ymax,
              zmin => $zmin, zmax => $zmax,
              zoom_xmin => $xmin, zoom_xmax => $xmax,
              zoom_ymin => $ymin, zoom_ymax => $ymax,
              zoom_zmin => $zmin, zoom_zmax => $zmax);
    return %instr;
}


sub transform {
    my ($comp, $x, $y, $z, $vx, $vy, $vz, $t, $ph1, $ph2) = @_;
    if(!$comp) {
        return ($x, $y, $z, $vx, $vy, $vz, $t, $ph1, $ph2);
    } else {
        my ($nx, $ny, $nz, $nvx, $nvy, $nvz, $nt, $nph1, $nph2);
        my @T = @{$transformations{$comp}};
        $x *= $magnification if $compdraw{$comp} && $compdraw{$comp}{'magX'};
        $y *= $magnification if $compdraw{$comp} && $compdraw{$comp}{'magY'};
        $z *= $magnification if $compdraw{$comp} && $compdraw{$comp}{'magZ'};
        $nx = $x*$T[3] + $y*$T[6] + $z*$T[9] + $T[0];
        $ny = $x*$T[4] + $y*$T[7] + $z*$T[10] + $T[1];
        $nz = $x*$T[5] + $y*$T[8] + $z*$T[11] + $T[2];
        $nvx = $vx*$T[3] + $vy*$T[6] + $vz*$T[9];
        $nvy = $vx*$T[4] + $vy*$T[7] + $vz*$T[10];
        $nvz = $vx*$T[5] + $vy*$T[8] + $vz*$T[11];
        return ($nx, $ny, $nz, $nvx, $nvy, $nvz, $t, $ph1, $ph2);
    }
}


sub get_inspect_pos {
    my ($inspect, @comps) = @_;
    return 0 unless $inspect;
    my $i;
    for($i = 0; $i < @comps; $i++) {
        return $i if $comps[$i] eq $inspect;
    }
    die "Error: Inspected component $inspect not part of instrument $sim ?";
}


sub read_neutron {
    my ($in) = @_;
    my (@x, @y, @z, @vx, @vy, @vz, @t, @ph1, @ph2, @ncomp);
    my ($st, $i);
    my ($comp, $numcomp, $EndFlag);

    $EndFlag = 0;
    $st = 0;
    $i = 0;
    $numcomp = 0;
    $EndFlag = 0;
    my $dropit = 1;                # Flag to drop uninteresting neutron states.
    while(<$in>) {
        if($st == 0 && /^ENTER:/) {
            # Neutron enters instrument.
            $st = 1;
        } elsif($st == 0 && /^STATE:/) {
            # State after leaving - should probably be removed in McStas.
            next;
        } elsif($st == 1 && /^COMP:\s*"([a-zA-Z0-9_]+)"\s*$/) {
            # Neutron enters component local coordinate system.
            $comp = $1;
            $numcomp++;
            $dropit = 1;        # Drop the first state (entry point).
        } elsif($st == 1 && (/^STATE:(.*)$/ || /^SCATTER:(.*)$/)) {
            # Neutron state.
            $dropit = 0, next if $dropit; # Skip entry point
            ($x[$i], $y[$i], $z[$i],
             $vx[$i], $vy[$i], $vz[$i],
             $t[$i], $ph1[$i], $ph2[$i]) = split ",", $1;
            ($x[$i], $y[$i], $z[$i],
             $vx[$i], $vy[$i], $vz[$i],
             $t[$i], $ph1[$i], $ph2[$i]) =
                 transform($comp, $x[$i], $y[$i], $z[$i],
                           $vx[$i], $vy[$i], $vz[$i],
                           $t[$i], $ph1[$i], $ph2[$i]);
            $ncomp[$i] = $comp;
            $i++;
        } elsif($st == 1 && /^ABSORB:/) {
            # Neutron was absorbed.
            next;                # No special action needed.
        } elsif($st == 1 && /^LEAVE:/) {
            # Neutron leaves instrument.
            $st = 2;
            last;
        } elsif (/^Detector:/){
          if (! $MCSTAS::mcstas_config{'PLOTTER'} =~ /scriptfile/i) {
            # Should only be done if finished, and not called with --save flag...
            # Also, can only be done if tcl/tl available
            if ($MCSTAS::mcstas_config{'TCLTK'} ne "no" && $EndFlag == 0) {
              my $main = new MainWindow;
              $main->Label(-text => "Simulation $sim ended."
                          )->pack;
              $main->Label(-text => 'Press Ok to terminate backend'
                          )->pack;
              $main->Button(-text => 'Ok',
                            -command => sub{destroy $main}
                           )->pack;
              MainLoop;
              $EndFlag = 1;
            }
          }
          print;
        } else {
            # Pass on any output not meant for us.
            print;
        }
    }
    exit unless $st == 2;        # Stop when EOF seen before neutron data end.

    my %neutron = ('x' => \@x, 'y' => \@y, z => \@z,
                   vx => \@vx, vy => \@vy, vz => \@vz,
                   t => \@t, ph1 => \@ph1, ph2 => \@ph2,
                   comp => \@ncomp, numcomp => $numcomp, EndFlag => $EndFlag);
    return %neutron
}


sub plot_components { # PGPLOT stuff only
    my ($rx, $ry, $rori, $rdis, $axis1, $axis2) = @_;
    my (@x, @y, @ori);
    my ($i, $col);

    @x = @$rx;
    @y = @$ry;
    @ori = @$rori;

    PGPLOT::pgsci(2);
    PGPLOT::pgline($#x + 1, \@x, \@y);
    PGPLOT::pgpt($#x + 1, \@x, \@y, 20);
    $col = 4;
    for($i = 0; $i <= $#components; $i++) {
        my $comp = $components[$i];
        PGPLOT::pgsci($col++);
        $col = 4 if $col > 15;
        if($compdraw{$comp}) {
            foreach $elem (@{$compdraw{$comp}{'elems'}}) {
                if($elem->{'type'} eq 'multiline') {
                    PGPLOT::pgline($elem->{'count'}, $elem->{$axis1}, $elem->{$axis2});
                }
            }
        } else {
            PGPLOT::pgsch(1.4);
            PGPLOT::pgpt(1, $x[$i], $y[$i], 26);
        }
    }
}


sub plot_neutron {
    my ($x, $y, $z, $vx, $vy, $vz, $comp) = @_;
    my ($i, $col, $oldcomp, $retval);
    if ($MCSTAS::mcstas_config{'PLOTTER'} =~ /McStas|PGPLOT/i) {
      # PGPLOT only
      PGPLOT::pgsci(3);
      PGPLOT::pgline(scalar(@$x), $x, $y);
      # Show component entry/exit points in same colour as respective component.
      # This should also be done w/ Matlab/Scilab...
      $i = 0;
      $col = 4;
      while($i < scalar(@$x)) {
        if(!defined($oldcomp) || $oldcomp cmp $comp->[$i]) {
          $oldcomp = $comp->[$i];
            PGPLOT::pgsci($col++);
          $col = 4 if $col > 15;
        }
        # Exit point.
        PGPLOT::pgpt(1, $x->[$i], $y->[$i], 17);
        $i++;
      }
    } elsif ($MCSTAS::mcstas_config{'PLOTTER'} =~ /Matlab/i) {
      # Matlab
      $retval=write_process("mcdisplay('Timeout');\n");
      $retval=write_process("mcdisplay('PlotNeutron',[@$x],[@$y],[@$z]);\n");
      return $retval;
    } elsif ($MCSTAS::mcstas_config{'PLOTTER'} =~ /Scilab/i) {
      # Scilab
      $retval=write_process("x=[];\n");
      $retval=write_process("y=[];\n");
      $retval=write_process("z=[];\n");
      $i=0;
      while($i < scalar(@$x)) {
        $tmp=$x->[$i];
        $retval=write_process("x=[x $tmp];\n");
        $tmp=$y->[$i];
        $retval=write_process("y=[y $tmp];\n");
        $tmp=$z->[$i];
        $retval=write_process("z=[z $tmp];\n");
        $i++;
      }
      $retval=write_process("PlotNeutron(x,y,z);\n");
      return $retval;
    }
}


sub show_comp_names { # PGPLOT stuff only
    my ($rinstr) = @_;
    my %instr = %$rinstr;
    my @comps = @{$instr{'comp'}};
    my $count = @comps;
    $count = 8 if $count < 8;
    my $col = 4;
    PGPLOT::pgsch(25/$count);
    my $i;
    for($i = 0; $i < @comps; $i++) {
        PGPLOT::pgsci($col++);
        $col = 4 if $col > 15;
        PGPLOT::pgmtxt('RV', 0.2, 1 - ($i+0.5)/$count, 0.0, $comps[$i]);
    }
}


sub reset_zoom {
    my ($rinstr, $vps) = @_;
    $rinstr->{'zoom_xmin'} = $rinstr->{'xmin'};
    $rinstr->{'zoom_xmax'} = $rinstr->{'xmax'};
    $rinstr->{'zoom_ymin'} = $rinstr->{'ymin'};
    $rinstr->{'zoom_ymax'} = $rinstr->{'ymax'};
    $rinstr->{'zoom_zmin'} = $rinstr->{'zmin'};
    $rinstr->{'zoom_zmax'} = $rinstr->{'zmax'};
    $zooming = 0;
}


sub do_zoom {
    my ($rinstr, $vps, $cx, $cy, $cx1, $cy1) = @_;
    my ($tmp, $a1, $a2);
    $tmp = $cx, $cx = $cx1, $cx1 = $tmp if $cx > $cx1;
    $tmp = $cy, $cy = $cy1, $cy1 = $tmp if $cy > $cy1;

    if($cx == $cx1 || $cy == $cy1) {
        print STDERR "Warning: bad zoom area.\n";
        return;
    }
    if($multi_view) {
        # Convert from screen coordinates to world coordinates.
        # First find which of the three views was choosen.
        if($cx < 0 && $cy < 1) {
            $cx = $cx + 1;
            $cx1 = $cx1 + 1;
            ($a1,$a2) = ("z", "y");
        } elsif($cx < 0 && $cy >= 1) {
            $cx = $cx + 1;
            $cx1 = $cx1 + 1;
            $cy = $cy - 1;
            $cy1 = $cy1 - 1;
            ($a1,$a2) = ("z", "x");
        } elsif($cx >= 0 && $cy >= 1) {
            $cy = $cy - 1;
            $cy1 = $cy1 - 1;
            ($a1,$a2) = ("x", "y");
        } else {
            print STDERR "Warning: bad zoom area.\n";
            return;
        }
        my $idx = "$a1-$a2";
        my $vpx0 = $vps->{$idx}{'VP'}[0];
        my $vpdx = $vps->{$idx}{'VP'}[1] - $vpx0;
        my $wx0 = $vps->{$idx}{'W'}[0];
        my $wdx = $vps->{$idx}{'W'}[1] - $wx0;
        my $vpy0 = $vps->{$idx}{'VP'}[2];
        my $vpdy = $vps->{$idx}{'VP'}[3] - $vpy0;
        my $wy0 = $vps->{$idx}{'W'}[2];
        my $wdy = $vps->{$idx}{'W'}[3] - $wy0;
        $cx = ($cx-$vpx0)/$vpdx*$wdx+$wx0;
        $cx1 = ($cx1-$vpx0)/$vpdx*$wdx+$wx0;
        $cy = ($cy-$vpy0)/$vpdy*$wdy+$wy0;
        $cy1 = ($cy1-$vpy0)/$vpdy*$wdy+$wy0;
    } else {
        ($a1, $a2) = ("z","x");
    }
    $rinstr->{"zoom_${a1}min"} = $cx;
    $rinstr->{"zoom_${a1}max"} = $cx1;
    $rinstr->{"zoom_${a2}min"} = $cy;
    $rinstr->{"zoom_${a2}max"} = $cy1;
    $zooming = 1;
}

sub write_process {
  my ($command) = @_;
  # $pid == 0 covers file output, Scilab/Matlab + Matlab on Win32 (OLE)
  if (!$pid eq 0) {
    (kill 0, $pid) || print STDERR "$plotter process terminated - ending...\n";
    return 2;
  }
  if ($Config{'osname'} eq 'MSWin32' && $MCSTAS::mcstas_config{'PLOTTER'} =~ /Matlab/i) {
    $ML->Execute($command);
    my $err = Win32::OLE::LastError();
    if ($err) {
      print STDERR "Matlab terminated! - Exiting.\n";
      return 2;
    }
  } else {
    # Simply write data to pipe/file
    print WRITER $command;
    return 1;
  }
}


sub plot_instrument {
    my ($noninteractive, $rinstr, $rneutron) = @_;
    my %instr = %$rinstr;
    my %neutron = %$rneutron;
    my $retval;
    # The following only relevant in PGPLOT mode
    if ($MCSTAS::mcstas_config{'PLOTTER'} =~ /McStas|PGPLOT/i) {
      my ($xmin, $xmax, $ymin, $ymax, $zmin, $zmax) =
        ($instr{'zoom_xmin'}, $instr{'zoom_xmax'}, $instr{'zoom_ymin'},
         $instr{'zoom_ymax'}, $instr{'zoom_zmin'}, $instr{'zoom_zmax'});
      my %vps;                        # Viewport/window setup.
      my ($vpx1,$vpx2,$vpy1,$vpy2,$wx1,$wx2,$wy1,$wy2);

      # PGPLOT::pgpage;        # start new page/panel
      PGPLOT::pgbbuf;        # begin buffer batch output

      # First show instrument from "above" (view in direction of y axis).

      PGPLOT::pgsci(1);
      PGPLOT::pgsch(1.4);
      PGPLOT::pgenv($zmin, $zmax, $xmin, $xmax, ($zooming ? 0 : 1), 0);
      PGPLOT::pglab("Z Axis [m]", "X Axis [m]", ($multi_view ? "Z-X view" : "Z-X view: $sim_cmd"));
      show_comp_names($rinstr);
      PGPLOT::pgsch(1.4);
      plot_components($instr{'z'}, $instr{'x'}, $instr{'ori'}, $instr{'dis'},
                      'Z', 'X');
      plot_neutron($neutron{'z'}, $neutron{'x'}, $neutron{'y'},
                   $neutron{'vz'}, $neutron{'vx'}, $neutron{'vy'},$neutron{'comp'});

      if($multi_view) {
        # Remember viewport setup for Z-X view.
        PGPLOT::pgqvp(0, $vpx1, $vpx2, $vpy1, $vpy2);
        PGPLOT::pgqwin($wx1, $wx2, $wy1, $wy2);
        $vps{'z-x'} = {VP => [$vpx1,$vpx2,$vpy1,$vpy2],
                       W => [$wx1,$wx2,$wy1,$wy2]};

        # Now show instrument viewed in direction of z axis.
        PGPLOT::pgsci(1);
        PGPLOT::pgsch(1.4);
        PGPLOT::pgenv($xmin, $xmax, $ymin, $ymax, ($zooming ? 0 : 1), 0);
        PGPLOT::pglab("X Axis [m]", "Y Axis [m]", "X-Y view");
        plot_components($instr{'x'}, $instr{'y'}, $instr{'ori'}, $instr{'dis'},
                        'X', 'Y');
        plot_neutron($neutron{'x'}, $neutron{'y'}, $neutron{'z'},
                     $neutron{'vx'}, $neutron{'vy'}, $neutron{'vz'} ,$neutron{'comp'});
        # Remember viewport setup for Z-X view.
        PGPLOT::pgqvp(0, $vpx1, $vpx2, $vpy1, $vpy2);
        PGPLOT::pgqwin($wx1, $wx2, $wy1, $wy2);
        $vps{'x-y'} = {VP => [$vpx1,$vpx2,$vpy1,$vpy2],
                       W => [$wx1,$wx2,$wy1,$wy2]};

        # Now show instrument viewed in direction of x axis.
        PGPLOT::pgsci(1);
        PGPLOT::pgsch(1.4);
        PGPLOT::pgenv($zmin, $zmax, $ymin, $ymax, ($zooming ? 0 : 1), 0);
        PGPLOT::pglab("Z Axis [m]", "Y Axis [m]", "Z-Y view");
        plot_components($instr{'z'}, $instr{'y'}, $instr{'ori'}, $instr{'dis'},
                        'Z', 'Y');
        plot_neutron($neutron{'z'}, $neutron{'y'}, $neutron{'x'},
                     $neutron{'vz'}, $neutron{'vy'}, $neutron{'vx'}, $neutron{'comp'});
        # Remember viewport setup for Z-Y view.
        PGPLOT::pgqvp(0, $vpx1, $vpx2, $vpy1, $vpy2);
        PGPLOT::pgqwin($wx1, $wx2, $wy1, $wy2);
        $vps{'z-y'} = {VP => [$vpx1,$vpx2,$vpy1,$vpy2],
                       W => [$wx1,$wx2,$wy1,$wy2]};

        # Set up viewport & window for mouse zoom.
        if ($multi_view) {
                PGPLOT::pgpanl(2,2);
                PGPLOT::pgsci(1);
                my $time=gmtime;
                PGPLOT::pgmtxt("t",0-1*1.2,0.0,0.0,"Date: $time");
                PGPLOT::pgmtxt("t",-2-1*1.2,0,0.0,"Simulation: ");
                PGPLOT::pgmtxt("t",-3-1*1.2,0.05,0.0,"$sim_cmd");
        }        # go to last panel
        else { PGPLOT::pgpanl(1,1); }                    # to activate full page
        PGPLOT::pgsvp(0,1,0,1);        # zoom for full page
        PGPLOT::pgswin(0,1,0,1);
      }
      PGPLOT::pgebuf;        # end buffer batch output

      return 0 if $noninteractive;

      # Now wait for a keypress in the graphics window.
      my ($cx, $cy, $cc);
      $cx = $cy = 0;
      PGPLOT::pgband(0, 0, 0, 0, $cx, $cy, $cc);
      if($cc =~ /[qQ]/) {
        return 2;                # Finished.
      } elsif($cc =~ /[pP]/) {        # B&W hardcopy.
        return 3;
      } elsif($cc =~ /[cC]/) {        # color hardcopy.
        return 4;
      } elsif($cc =~ /[gG]/) {        # GIF hardcopy.
        return 5;
      } elsif($cc =~ /[nN]/) {        # PNG hardcopy.
        return 6;
     } elsif($cc =~ /[hH]/) {        # Help
        print STDERR "McDisplay: q=quit, h=help\n";
        print STDERR "    output p=ps,   c=color ps, g=gif n=png\n";
        print STDERR "      zoom x=reset,z or d=zoom\n";
      } elsif($cc =~ /[zZdD]/) {        # Zoom.
        my ($cx1, $cy1, $cc1) = (0, 0, 0);
        PGPLOT::pgband(2,0,$cx,$cy,$cx1,$cy1,$cc1);
        do_zoom($rinstr, \%vps, $cx, $cy, $cx1, $cy1);
        return 1;
      } elsif($cc =~ /[xX]/) {        # Reset zoom.
              PGPLOT::pgpanl(1,1);
        PGPLOT::pgperas;
        if ($multi_view) {
                PGPLOT::pgpanl(1,2);
                PGPLOT::pgperas;
                PGPLOT::pgpanl(2,1);
                PGPLOT::pgperas;
                PGPLOT::pgpanl(2,2);
                PGPLOT::pgperas;
        }
        reset_zoom($rinstr, \%vps);
        return 1;
      }
    } else {
      # Leave further checks for plot_neutron
      $retval=plot_neutron($neutron{'z'}, $neutron{'x'}, $neutron{'y'},
                   $neutron{'vz'}, $neutron{'vx'}, $neutron{'vy'}, $neutron{'comp'});
      if ($retval==2) {
        return $retval;
      }
    }
    return 0;                        # Default: do not repeat this neutron.
}

sub get_device { # PGPLOT stuff only
    my ($what) = @_;
    my $dev;

    if (defined(&dev)) { $dev = dev($what); }
    else { $dev = PGPLOT::pgopen($what); }
    return $dev if $dev < 0;
    if($multi_view) {
    # We use a 2x2 display format to view the instrument from three angles.
        PGPLOT::pgsubp(2, 2);
    } else {
        # We use a 1x1 display for detail.
        PGPLOT::pgsubp(1, 1);
    }
    return $dev;
}


# Test the code.

# Check command line arguments.

undef $inspect;
undef $first;
undef $last;
undef $save;
undef $direct_output;
undef $sim_cmd;
undef $sim;
my $plotter;
undef $file_output;
my $int_mode=0; # interactive mode(0), non interactive (1)
my $i;
my $start_scilab=0;
my $show_help=0;

$plotter = $MCSTAS::mcstas_config{'PLOTTER'};

for($i = 0; $i < @ARGV; $i++) {
    if(($ARGV[$i] eq "-m") || ($ARGV[$i] eq "--multi")) {
        $multi_view = 1;
    } elsif($ARGV[$i] =~ /--help|-h$/) {
        $show_help=1;
    } elsif(($ARGV[$i] =~ /^-z([-0-9+.eE]+)$/) ||
            ($ARGV[$i] =~ /^--zoom=([-0-9+.eE]+)$/)) {
        $magnification = ($1 == 0 ? 1 : $1);
    } elsif(($ARGV[$i] eq "-gif") || ($ARGV[$i] eq "-ps") ||
            ($ARGV[$i] eq "-fig") || ($ARGV[$i] eq "-scg") ||
            ($ARGV[$i] eq "-psc") || ($ARGV[$i] eq "-png") || ($ARGV[$i] eq "-ppm")) {
        $direct_output = $ARGV[$i];
        $save = 1;
        $int_mode = 1;
    } elsif(($ARGV[$i] =~ /^-i([a-zA-Z0-9_]+)$/) ||
            ($ARGV[$i] =~ /^--inspect=([a-zA-Z0-9_]+)$/)) {
        $inspect = $1;
    } elsif($ARGV[$i] =~ /^--first=([a-zA-Z0-9_]+)$/) {
        $first = $1;
    } elsif($ARGV[$i] =~ /^--last=([a-zA-Z0-9_]+)$/) {
        $last = $1;
    } elsif($ARGV[$i] eq "--save") {
        $save = 1;
    } elsif(($ARGV[$i] =~ /^-p([a-zA-Z0-9_]+)$/) ||
              ($ARGV[$i] =~ /^--plotter=([a-zA-Z0-9_\"]+)$/) ||
              ($ARGV[$i] =~ /^--format=([a-zA-Z0-9_\"]+)$/)) {
        $plotter = $1;
   } elsif(($ARGV[$i] =~ /^-f([a-zA-Z0-9_\-\/\ \.\:\"]+)$/) ||
              ($ARGV[$i] =~ /^--file=([a-zA-Z0-9_\-\/\ \.\:]+)$/)) {
        $file_output = $1;
   } else {
        if (defined($sim_cmd)) { push @cmdline, $ARGV[$i]; }
        else {
          $sim_cmd = $ARGV[$i];
          $sim=$sim_cmd;
          # Remove trailing .out or .exe extension
          $sim=~ s|.out\Z||;
          $sim=~ s|.exe\Z||;
          $sim=~ s|.instr\Z||;
        }
   }
}
if ($show_help) { undef $sim_cmd; }
die "Usage: mcdisplay [-mzipfh][-gif|-ps|-psc] Instr.out [instr_options] params
 -h        --help            Show this help
 -m        --multi           Show the three instrument side views
 -zZF      --zoom=ZF         Show zoomed view by factor ZF
 -iCOMP    --inspect=COMP    Show only trajectories reaching component COMP
 -pPLOTTER --plotter=PLOTTER Output graphics using {PGPLOT,Scilab,Matlab}
 -fFNAME   --file=FNAME      Output graphics commands to file FNAME
                             (Only used when PLOTTER = {Scilab, Matlab})
           --first=COMP      First component to visualize {Scilab, Matlab}
           --last=COMP       Last component to visualize {Scilab, Matlab}
           --save            Output a Scilab/Matlab figure file and exit
                             (Filename is Instr.scf / Instr.fig). Figure
                             files are used by mcgui.pl for visualising the
                             instrument. With PGPLOT, --save is nonfunctional.
 -gif|-ps|-psc               Export figure as gif/b&w ps/color ps and exit
 When using -ps -psc -gif, the program writes the hardcopy file and exits.
 SEE ALSO: mcstas, mcdoc, mcplot, mcrun, mcgui, mcresplot, mcstas2vitess
 DOC:      Please visit http://www.mcstas.org/\n"
 unless $sim_cmd;

if ($sim_cmd =~ m'\.instr$') # recompile .instr if needed
{ my @ccopts=();
  $sim_cmd = get_out_file($sim_cmd, 0, @ccopts); }


# Check value of $plotter and $file_output variables, set
# $MCSTAS::mcstas_config{'PLOTTER'} with scriptfile keyword
if ($file_output) { $plotter .= "_scriptfile"; }

if ($plotter =~ /Scilab/i && not $plotter =~ /scriptfile/i) {
  # Check for Win32, only file_output possible :(
  if ($Config{'osname'} eq MSWin32) {
    $file_output="mcdisplay_commands.sci";
    print STDERR "\n******************************************************\n";
    print STDERR "Sorry, Scilab only possible using file output on Win32\n\n";
    print STDERR "Outputting to file $file_output\n";
    print STDERR "******************************************************\n\n";
    print STDERR "When done, I will start a scilab for you to view the file...\n";
    $start_scilab=1;
    $plotter="Scilab_scriptfile";
  }
}

if ($plotter =~ /scriptfile/i && not $file_output) {
  $file_output="mcdisplay_commands";
  if ($plotter =~ /Matlab/i) { $file_output .=".m"; }
  elsif ($plotter =~ /Scilab/i) { $file_output .=".sci"; }
  print STDERR "Outputting to file $file_output\n";
}

# Final PLOTTER check, is PGPLOT wanted but not possible?
# - Ask user to rerun / set other default
if ($plotter =~ /McStas|PGPLOT/i && $MCSTAS::mcstas_config{'PGPLOT'} eq "no") {
  print STDERR "\n******************************************************\n";
  print STDERR "Default / selected PLOTTER is PGPLOT - Problems:\n\n";
  print STDERR "PGPLOT.pm not found on Perl \@INC path\n\nSolution:\n\n";
  print STDERR "1) Install pgplot + pgperl packages (Unix/Linux/Cygwin) \n";
  print STDERR "2) Rerun mcdisplay with -p/--plotter set to Scilab/Matlab \n";
  print STDERR "3) Modify $MCSTAS::perl_dir/mcstas_config.perl\n";
  print STDERR "   to set a different default plotter\n\n";
  print STDERR "******************************************************\n\n";
  die "mcdisplay: PGPLOT problems...\n";
}

$MCSTAS::mcstas_config{'PLOTTER'} = $plotter;

my $pg_devname = "/xserv";
# Only set up PGPLOT stuff if needed
if ($plotter =~ /McStas|PGPLOT/i) { # PGPLOT is plotter!
  if ($int_mode == 1)
    {
      my $ext  = "ps";
      my $type = "ps";
      if($direct_output eq "-gif") { $ext="gif"; $type="gif"; }
      elsif($direct_output eq "-png") { $ext="png"; $type="png"; }
      elsif($direct_output eq "-psc") { $type="cps"; }
      $pg_devname = "$sim_cmd.$ext/$type";
    }
  $global_device = get_device($pg_devname);
  if($global_device < 0) {
    print STDERR "mcdisplay: Failed to open PGPLOT device $pg_devname\n";
    exit 1;
  }
  my $seq = "";                        # Sequence number for multiple hardcopy
  PGPLOT::pgask(0);
} elsif ($plotter =~ /Matlab/i && $plotter =~ /scriptfile/i) {
  # Matlab w/FILE is plotter - open a file handle
  open(WRITER, "> $file_output");
  $pid=0;
} elsif ($plotter =~ /Matlab/i) {
  # Matlab is plotter - open a pipe / OLE connection
  if ($Config{'osname'} eq 'MSWin32') {
    $pid=0;
    $ML = Win32::OLE->new('Matlab.Application') || die "mcdisplay: Could not start Matlab\n";
  } else {
    my $cmd = "$MCSTAS::mcstas_config{'MATLAB'} $MCSTAS::mcstas_config{'MATLAB_COMMAND'} > /dev/null";
    $pid=open2(READER,WRITER, $cmd) || die "mcdisplay: Could not start Matlab\n$cmd\n";
  }
  print STDERR "Opened up pipe to matlab - pid is $pid\n";
  print STDERR "Building Matlab INSTRUMENT struct, please have patience...\n";
} elsif ($plotter =~ /Scilab/i && $plotter =~ /scriptfile/i) {
  # Scilab w/FILE is plotter - open a file handle
  open(WRITER, "> $file_output");
  $pid=0;
} elsif ($plotter =~ /Scilab/i) {
  # Scilab is plotter - open a pipe
  my $cmd = "$MCSTAS::mcstas_config{'SCILAB'} $MCSTAS::mcstas_config{'SCILAB_COMMAND'} > /dev/null";
  $pid=open2(READER,WRITER, $cmd) || die "Could not start Scilab\n$cmd\n";
  print STDERR "Opened up pipe to scilab - pid is $pid\n";
  print STDERR "Building Scilab INSTRUMENT struct, please have patience...\n";
}

my ($numcomp, %neutron, %instr);

$args = join(" ", @cmdline);
$cmdline = "$sim_cmd --trace --no-output-files $args";
printf STDERR "Starting simulation '$cmdline' ...\n";
open(IN, "$cmdline |") || die "mcdisplay: Could not run simulation\n";

$numcomp = read_instrument(IN);
$inspect_pos = get_inspect_pos($inspect, @components);
%instr = make_instrument;

if ($int_mode == 0 && $plotter =~ /McStas|PGPLOT/i)  { printf STDERR "McDisplay/PGPLOT: Press H key for help.\n"; }
if ($plotter =~ /Matlab/i && not $plotter =~ /scriptfile/i) { print STDERR "Matlab INSTRUMENT done, starting gui....\n"; }
if ($plotter =~ /Scilab/i && not $plotter =~ /scriptfile/i) { print STDERR "Scilab INSTRUMENT done, starting gui....\n"; }

while(!eof(IN)) {
    %neutron = read_neutron(IN);
    if ($start_scilab == 1) {
    # This only happens on Win32 (runscilab.exe), and we know the filename too...
        my $runscilab = "$MCSTAS::mcstas_config{'SCILAB'}";
        my $pid = open(SCILAB,"$runscilab -f mcdisplay_commands.sci|");
        while(!eof(SCILAB)) {
            # Do nothing...
        }
        $start_scilab = 0;
    }
    next if $neutron{'numcomp'} <= $inspect_pos;

    my $ret;
    do {
        $ret = plot_instrument($int_mode, \%instr, \%neutron);
        if ($plotter =~ /McStas|PGPLOT/i) {
          if ($int_mode == 1) { $ret =2; print STDERR "Wrote \"$pg_devname\"\n"; }
          if($ret == 3 || $ret == 4 || $ret == 5 || $ret == 6) {
            my $ext="ps";
            my $type = $ret == 3 ? "ps" : "cps";
            if($ret == 5) { $type = "gif"; $ext="gif"; }
            if($ret == 6) { $type = "png"; $ext="png"; }
            my $tmp_pg_devname = "$sim_cmd$seq.$ext/$type";
            my $tmpdev=0;
            $tmpdev = get_device($tmp_pg_devname);
            if($tmpdev < 0) {
              print STDERR "mcdisplay: Warning: could not open PGPLOT output \"$tmp_pg_devname\" for hardcopy output\n";
            } else {
              plot_instrument(1, \%instr, \%neutron);
              print STDERR "Wrote \"$tmp_pg_devname\"\n";
              ++$seq;
	    }
            if (defined(&close_window)) { 
		close_window($tmpdev); 
	    } else {
		PGPLOT::pgclos();
	    }
	    PGPLOT::pgslct($global_device);
#	    PGPLOT::pgask(0);
          }
        }
    } while($ret != 0 && $ret != 2);
    last if $ret == 2;
  }
close(IN);

# Properly close any open files etc.
if ($plotter =~ /McStas|PGPLOT/i) {
  if (defined(&close_window)) { close_window(); }
  else {  PGPLOT::pgclos(); }
} else {
  close(WRITER);
}
