root/branches/mcstas-1.x/mcdisplay.pl

Revision 2699, 48.5 KB (checked in by pkwi, 2 years ago)

Added files from mcstas-1.12a
  • Property svn:executable set to *
Line 
1#! /usr/bin/perl
2# Removed -w to get rid of "used only once" warnings
3
4# In emacs, please make this -*- perl -*- mode. Thanks.
5
6# Enhanced mcdisplay script for --trace output from McStas simulation.
7#
8# Implements graphic display using either
9# PGPLOT
10# Scilab
11# Matlab
12#
13# PW 20030320
14#
15#   This file is part of the McStas neutron ray-trace simulation package
16#   Copyright (C) 1997-2004, All rights reserved
17#   Risoe National Laborartory, Roskilde, Denmark
18#   Institut Laue Langevin, Grenoble, France
19#
20#   This program is free software; you can redistribute it and/or modify
21#   it under the terms of the GNU General Public License as published by
22#   the Free Software Foundation; version 2 of the License.
23#
24#   This program is distributed in the hope that it will be useful,
25#   but WITHOUT ANY WARRANTY; without even the implied warranty of
26#   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
27#   GNU General Public License for more details.
28#
29#   You should have received a copy of the GNU General Public License
30#   along with this program; if not, write to the Free Software
31#   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
32
33
34
35
36
37
38# Config module needed for Win32 vs unix setup.
39# PW 20030320
40use Config;
41
42BEGIN {
43  # default configuration (for all high level perl scripts)
44  if($ENV{"MCSTAS"}) {
45    $MCSTAS::sys_dir = $ENV{"MCSTAS"};
46  } else {
47    if ($Config{'osname'} eq 'MSWin32') {
48      $MCSTAS::sys_dir = "c:\\mcstas\\lib";
49    } else {
50      $MCSTAS::sys_dir = "/usr/local/lib/mcstas";
51    }
52  }
53  $MCSTAS::perl_dir = "$MCSTAS::sys_dir/tools/perl";
54  $MCSTAS::perl_modules = "$MCSTAS::perl_dir/modules";
55
56  # custom configuration (this script)
57
58  # If this is Win32, load OLE related modules -
59  # we can not talk to Matlab using pipe on Win32 :(
60  # PW 20030314
61  if ($Config{'osname'} eq 'MSWin32') {
62    require Win32::OLE;
63    require Win32::OLE::Variant;
64  }
65}
66
67use lib $MCSTAS::perl_dir;
68use lib $MCSTAS::perl_modules;
69require "mcstas_config.perl";
70
71# Overload with user's personal config
72if ($ENV{"HOME"} && -e $ENV{"HOME"}."/.mcstas/mcstas_config.perl") {
73  require $ENV{"HOME"}."/.mcstas/mcstas_config.perl";
74}
75
76require "mcrunlib.pl";
77# IPC can probably be used safely, exists on sysv type systems,
78# linux, Win32. Will investigate further regarding portability.
79# PW 20030404
80use IPC::Open2;
81use Math::Trig;
82
83$magnification = 1;
84$zooming = 0;
85
86my (%transformations, @components);
87
88
89sub read_instrument {
90    my ($in) = @_;
91    my ($st, $comp);
92    my $compheader;
93    $st = 0;
94    @components = ();
95    while(<$in>) {
96        if($st == 0 && /^INSTRUMENT:/) {
97            # Start of instrument description.
98            $st = 1;
99            if ($MCSTAS::mcstas_config{'PLOTTER'} =~ /Matlab/i) {
100              # Initialize matlab struct...
101              write_process("addpath('$MCSTAS::sys_dir/tools/matlab');\n");
102              write_process("mcdisplay('Init');\n");
103              write_process("global INSTRUMENT;\n");
104              write_process("INSTRUMENT.descr='$sim';\n");
105              # Possibly, set firstcomp + lastcomp
106              if ($first) { write_process("INSTRUMENT.firstcomp='$first';\n"); }
107              if ($last)  { write_process("INSTRUMENT.lastcomp ='$last';\n");  }
108            }
109            if ($MCSTAS::mcstas_config{'PLOTTER'} =~ /Scilab/i) {
110              # Initialize scilab struct...
111              write_process("exec('$MCSTAS::sys_dir/tools/scilab/mcdisplay.sci',-1);\n");
112              write_process("INSTRUMENT.descr='$sim';\n");
113              # Possibly, set firstcomp + lastcomp
114              if ($first) { write_process("INSTRUMENT.firstcomp='$first';\n"); }
115              if ($last)  { write_process("INSTRUMENT.lastcomp ='$last';\n");  }
116              if ($save) {
117                if ($direct_output) { write_process("INSTRUMENT.save_format='$direct_output';\n"); }
118                else  { write_process("INSTRUMENT.save_format='-scg';\n"); }
119                write_process("INSTRUMENT.save=1;\n");
120              } else { write_process("INSTRUMENT.save=0;\n"); }
121            }
122            if ($MCSTAS::mcstas_config{'PLOTTER'} =~ /VRML/i) {
123                # Default viewpoint, 10 meters along z.
124                write_process("#VRML V2.0 utf8\n
125# Format: VRML 2.0\n
126# Output from mcdisplay from the McStas package, see http://www.mcstas.org
127# use freeWRL, openvrml, vrmlview, CosmoPlayer, Cortona, Octaga... to view file
128#\n# Instrument used was $sim_cmd. Full cmdline was:\n#
129# mcdisplay @ARGV
130#\n# Please rerun instrument with -i option to get more info.\n#\n
131WorldInfo {
132  title \"McStas: $sim_cmd instrument\"
133  info [ \"URL:    http://www.mcstas.org/\"
134    \"Editor: mcdisplay @ARGV\"
135    \"Creator:$sim_cmd simulation (McStas)\" ]
136}
137Viewpoint {
138  description \"Default\"
139   position 0 0.2 -1
140  orientation 0 1 0 3.14
141  jump FALSE
142}
143Background {
144  skyAngle [ 1.57 1.57]
145  skyColor [0 0 1, .1 .1 .1, 0.1 0 0]
146}\n");
147                # Definition of Origin + coordinate system arrows
148                write_process("# Sphere at the origin
149Shape {
150  appearance Appearance {
151  material Material {
152  diffuseColor 1.0 1.0 0.0
153  transparency 0.5 } }
154  geometry Sphere { radius 0.01 }
155}
156# Axis-parallel arrows of length 1 metre
157DEF ARROW Group {
158children [
159Transform {
160  translation 0 0.5 0
161  children [
162  Shape {
163  appearance DEF ARROW_APPEARANCE Appearance {
164  material Material {
165  diffuseColor .3 .3 1
166  emissiveColor .1 .1 .33
167  }
168  }
169  geometry Cylinder {
170  bottom FALSE
171  radius .005
172  height 1
173  top FALSE
174  } } ] }
175Transform {
176  translation 0 1 0
177  children [
178  DEF ARROW_POINTER Shape {
179  geometry Cone {
180  bottomRadius .05
181  height .1
182  }
183  appearance USE ARROW_APPEARANCE
184  } ] }
185] }
186# the arrow along X axis
187Transform {
188  translation 0 0 0
189  rotation 1 0 0 1.57
190  children [
191  Group {
192  children [
193  USE ARROW
194  ] } ] }
195# the arrow along Z axis
196Transform {
197  translation 0 0 0
198  rotation 0 0 1 -1.57
199  children [
200  Group {
201  children [
202  USE ARROW
203  ] } ] }
204# the Y label (which is vertical)
205DEF Y_Label Group {
206children [
207Transform {
208  translation 0 1 0
209  children [
210  Billboard {
211  children [
212  Shape {
213  appearance DEF LABEL_APPEARANCE Appearance {
214  material Material {
215  diffuseColor 1 1 .3
216  emissiveColor .33 .33 .1
217  } }
218  geometry Text {
219  string [\"y\" ]
220  fontStyle FontStyle {  size .2 }
221  } } ] } ] }
222] }
223# the X label
224DEF X_Label Group {
225children [
226Transform {
227  translation 1 0 0
228  children [
229  Billboard {
230  children [
231  Shape {
232  appearance DEF LABEL_APPEARANCE Appearance {
233  material Material {
234  diffuseColor 1 1 .3
235  emissiveColor .33 .33 .1
236  } }
237  geometry Text {
238  string [\"x\"]
239  fontStyle FontStyle {  size .2 }
240  } } ] } ] }
241] }
242# the Z label
243DEF Z_Label Group {
244children [
245Transform {
246  translation 0 0.2 1
247  children [
248  Billboard {
249  children [
250  Shape {
251  appearance DEF LABEL_APPEARANCE Appearance {
252  material Material {
253  diffuseColor 1 1 .3
254  emissiveColor .33 .33 .1
255  } }
256  geometry Text {
257  string [\"z\"]
258  fontStyle FontStyle {  size .2 }
259  } } ] } ] }
260] }
261# The text information (header data )
262DEF Header Group {
263children [
264Transform {
265  translation 0 1.2 0
266  children [
267  Billboard {
268  children [
269  Shape {
270  appearance Appearance {
271  material Material {
272  diffuseColor .9 0 0
273  emissiveColor .9 0 0 }
274  }
275  geometry Text {
276  string [ \"McStas: $sim_cmd\" ]
277  fontStyle FontStyle {
278  style \"BOLD\"
279  size .2
280  } } } ] } ] }
281] }
282\n### Instrument begins here: ###\n");
283
284            }
285        } elsif($st == 1 && /^COMPONENT:\s*"([a-zA-Z0-9_]+)"\s*/) {
286            $comp = $1;
287            @components = (@components, $comp);
288            if ($MCSTAS::mcstas_config{'PLOTTER'} =~ /Matlab/i) {
289              # Initialize components in matlab struct:
290              write_process("INSTRUMENT.name{length(INSTRUMENT.name)+1}='$comp';\n");
291            }
292            if ($MCSTAS::mcstas_config{'PLOTTER'} =~ /Scilab/i) {
293              # Initialize components in scilab struct:
294              write_process("setcomponent('$comp');\n");
295            }
296            if ($MCSTAS::mcstas_config{'PLOTTER'} =~ /vrml/i) {
297                $compheaders{$comp}="\n########################".
298                    " $comp ".
299                    "########################";
300            }
301        } elsif($st == 1 && /^POS:(.*)$/) {
302            my @T;
303            @T = split ",", $1;
304            $transformations{$comp} = \@T;
305            $compcnt=scalar(@components);
306            if ($MCSTAS::mcstas_config{'PLOTTER'} =~ /Matlab/i) {
307              # Component position for matlab struct:
308              write_process("INSTRUMENT.$comp.T=ReshapeTransform([@T]);\n");
309              write_process("INSTRUMENT.$comp.j=$compcnt;\n");
310              write_process("INSTRUMENT.$comp.K=cell(0);\n");
311            }
312            if ($MCSTAS::mcstas_config{'PLOTTER'} =~ /Scilab/i) {
313              # Component position for scilab struct:
314              write_process("setposition([@T]);\n");
315            }
316            if ($MCSTAS::mcstas_config{'PLOTTER'} =~ /vrml/i) {
317                if($T[0]!=0 or $T[1]!=0 or $T[2]!=0){$transforms{$comp}.= "translation $T[0] $T[1] $T[2]\n";}
318                my $angle = acos(($T[3]+$T[7]+$T[11]-1)/2);
319                my $d21=$T[8]-$T[10]; my $d02=$T[9]-$T[5]; my $d10=$T[4]-$T[6];
320                my $d=sqrt($d21*$d21+$d02*$d02+$d10*$d10);
321
322                if($d!=0){$transforms{$comp}.= "rotation ".$d21/$d.' '.$d02/$d.' '.$d10/$d." $angle\n";}
323                $nbcomp1++;
324                if($transforms{$comp})
325                {
326                    $compheaders{$comp}="$compheaders{$comp}".
327                        "\nTransform {\n".
328                        "$transforms{$comp}".
329                        "children [";
330                }
331                $compheaders{$comp}="$compheaders{$comp}".
332                    "\nShape {\nappearance Appearance {\n\t".
333                    "material Material { emissiveColor ";
334                my $color = vrml_setcolor($nbcomp2,$nbcomp1);
335                $compheaders{$comp}="$compheaders{$comp}".
336                    "$color".
337                    "}}\n".
338                    "geometry IndexedLineSet {\ncoord Coordinate {\npoint [\n";
339                $nbcomp2++;
340                $comp = "";
341            }
342
343        } elsif($st == 1 && /^MCDISPLAY: start$/) {
344            $st = 2;                # Start of component graphics representation
345        } elsif($st == 2 && /^MCDISPLAY: component ([a-zA-Z0-9_]+)/) {
346            $comp = $1;
347            $compdraw{$comp} = {};
348            $compdraw{$comp}{'elems'} = [];
349            #$compdraw{$comp}{'header'} = $compheader;
350            if ($MCSTAS::mcstas_config{'PLOTTER'} =~ /Scilab/i) {
351              # Initialize component variable etc. in scilab
352              write_process("trace_component('$comp');\n");
353            }
354        } elsif($st == 2 && /^MCDISPLAY: magnify\('([xyz]*)'\)$/) {
355            my $mag = $1;
356            $compdraw{$comp}{'magX'} = 1 if $mag =~ /x/i;
357            $compdraw{$comp}{'magY'} = 1 if $mag =~ /y/i;
358            $compdraw{$comp}{'magZ'} = 1 if $mag =~ /z/i;
359        } elsif($st == 2 && /^MCDISPLAY: multiline\(([0-9]+),([^()\n]+)\)$/) {
360            my $count = $1;
361            my @coords = split ',', $2;
362            push @{$compdraw{$comp}{'elems'}},
363                {type => 'multiline',
364                 count => $count,
365                 coords => \@coords};
366            if ($MCSTAS::mcstas_config{'PLOTTER'} =~ /Matlab/i) {
367              # Line elements for Matlab struct
368              write_process("coords=[@coords];\n");
369              write_process("x=coords(1:3:length(coords));\n");
370              write_process("y=coords(2:3:length(coords));\n");
371              write_process("z=coords(3:3:length(coords));\n");
372              write_process("coords=[x;y;z;1+0*z];\n");
373              write_process("INSTRUMENT.$comp.K{size(INSTRUMENT.$comp.K,2)+1}=coords;\n");
374            }
375            if ($MCSTAS::mcstas_config{'PLOTTER'} =~ /Scilab/i) {
376              # Line elements for Scilab struct
377              write_process("multiline([$1 @coords]);\n");
378            }
379        } elsif($st == 2 &&
380                /^MCDISPLAY: circle\('([xyzXYZ]{2})',([-+0-9.eE]+),([-+0-9.eE]+),([-+0-9.eE]+),([-+0-9.eE]+)\)$/) {
381            my ($plane,$x,$y,$z,$r) = ($1,$2,$3,$4,$5);
382            # Make a circle using a 25-order multiline.
383            my @coords = ();
384            my $i;
385            for($i = 0; $i <= 24; $i++) {
386                my $a = $r*cos(2*3.1415927/24*$i);
387                my $b = $r*sin(2*3.1415927/24*$i);
388                my ($x1,$y1,$z1) = ($x,$y,$z);
389                if($plane =~ /xy|yx/i) {
390                    $x1 += $a;
391                    $y1 += $b;
392                } elsif($plane =~ /xz|zx/i) {
393                    $x1 += $a;
394                    $z1 += $b;
395                } elsif($plane =~ /yz|zy/i) {
396                    $y1 += $a;
397                    $z1 += $b;
398                } else {
399                    die "mcdisplay: Bad plane specifier in circle: '$plane'";
400                }
401                push @coords, $x1, $y1, $z1;
402            }
403            push @{$compdraw{$comp}{'elems'}},
404                {type => 'multiline',
405                 count => 25,
406                 coords => \@coords};
407            if ($MCSTAS::mcstas_config{'PLOTTER'} =~ /Matlab/i) {
408              # Line elements for Matlab struct, circle representation
409              write_process("coords=[@coords];\n");
410              write_process("x=coords(1:3:length(coords));\n");
411              write_process("y=coords(2:3:length(coords));\n");
412              write_process("z=coords(3:3:length(coords));\n");
413              write_process("coords=[x;y;z;1+0*z];\n");
414              write_process("INSTRUMENT.$comp.K{size(INSTRUMENT.$comp.K,2)+1}=coords;\n");
415            }
416            if ($MCSTAS::mcstas_config{'PLOTTER'} =~ /Scilab/i) {
417              # Line elements for Scilab struct, circle representation
418              write_process("circle('$plane',$x,$y,$z,$r);\n");
419            }
420        } elsif($st == 2 && /^MCDISPLAY: end$/) {
421            $st = 1;  # End of component graphics representation
422            if ($MCSTAS::mcstas_config{'PLOTTER'} =~ /Matlab/i) {
423              # Matlab 'End of instrument'
424              write_process("mcdisplay('Load');\n");
425              write_process("PlotInstrument('init');\n");
426              # Check if we were called with --save option, output matlab figure if so...
427              if ($save) {
428                # Clone the graph to another window...
429                write_process("ax=gca;\n");
430                write_process("h=figure('numbertitle','off','name','$sim McStas Instrument')\n;");
431                write_process("copyobj(ax,h);\n");
432                write_process("saveas(h,'$sim.fig','fig');\n");
433                write_process("delete(h);\n");
434                write_process("delete(INSTRUMENT.fig);\n");
435                write_process("exit;\n");
436              } else {
437                write_process("set(INSTRUMENT.fig, 'closerequestfcn','exit;');\n");
438                write_process("wait(INSTRUMENT.fig);\n");
439              }
440            }
441            if ($MCSTAS::mcstas_config{'PLOTTER'} =~ /Scilab/i) {
442              # Scilab 'End of instrument'
443              write_process("endtrace();\n");
444            }
445            if ($MCSTAS::mcstas_config{'PLOTTER'} =~ /vrml/i) {
446                foreach $comp (@components) {
447                    write_process("$compheaders{$comp}");
448                    my @multilineSize =();
449                    foreach $elem (@{$compdraw{$comp}{'elems'}}) {
450                        push @multilineSize, $elem->{'count'};
451                        my @coords = @{$elem->{'coords'}};
452                        my $i=0;
453                        my $coordstring="";
454                        foreach $coord (@coords) {
455                            if ($i<2) {
456                                $coordstring="$coordstring$coord ";
457                                $i++;
458                            } else {
459                                $i=0;
460                                $coordstring="$coordstring$coord,";
461                            }
462                        }
463                        write_process("$coordstring\n");
464                    }
465                    write_process("]}\n");
466                    write_process("coordIndex [\n");
467
468                    my $c=0,$i,$j;
469                    my $m = join('/',@multilineSize);
470                    foreach $i (@multilineSize)
471                    {
472                        for($j=0; $j<$i; $j++)
473                        {
474                            write_process("$c,");
475                            $c++;
476                        }
477                        write_process("-1,\n");
478                    }
479                    write_process("]}}\n");
480                    if($transforms{$comp})
481                    {
482                        write_process("]}\n");
483                    }
484                    my @T=@{$transformations{$comp}};
485
486                    write_process("Viewpoint {\n".
487                                  "description \"$comp\"".
488                                  " position".$T[0].' '.$T[1].' '.$T[2]."\n".
489                                  "orientation 0 1 0  3.14\n".
490                                  "jump FALSE ".
491                                  "}\n");
492                }
493            }
494
495        } elsif($st == 1 && /^INSTRUMENT END:/) {
496            $st = 100;
497            last;
498        } else {
499            print;
500        }
501    }
502    exit if($st != 100);        # Stop when EOF seen before instrument end.
503    return $#components + 1;
504
505}
506
507
508sub make_instrument {
509    my (@x, @y, @z, @ori, @dis, @comp);
510    my ($i, $c, %instr);
511    my ($xmin, $xmax, $ymin, $ymax, $zmin, $zmax);
512
513    $i = 0;
514    foreach $c (@components) {
515        my (@T, @U);
516        @T = @{$transformations{$c}};
517        $x[$i] = $T[0];
518        $xmin = $x[$i] if !$xmin || $x[$i] < $xmin;
519        $xmax = $x[$i] if !$xmax || $x[$i] > $xmax;
520        $y[$i] = $T[1];
521        $ymin = $y[$i] if !$ymin || $y[$i] < $ymin;
522        $ymax = $y[$i] if !$ymax || $y[$i] > $ymax;
523        $z[$i] = $T[2];
524        $zmin = $z[$i] if !$zmin || $z[$i] < $zmin;
525        $zmax = $z[$i] if !$zmax || $z[$i] > $zmax;
526
527        @U = ($T[3], $T[4], $T[5], $T[6], $T[7], $T[8], $T[9], $T[10], $T[11]);
528        $ori[$i] = \@U;
529        $comp[$i] = $c;
530        # Now transform coordinates for component graphics representations.
531        if($compdraw{$c}) {
532            my $magX = $compdraw{$c}{'magX'};
533            my $magY = $compdraw{$c}{'magY'};
534            my $magZ = $compdraw{$c}{'magZ'};
535            foreach $elem (@{$compdraw{$c}{'elems'}}) {
536                if($elem->{'type'} eq 'multiline') {
537                    my @coords = @{$elem->{'coords'}};
538                    my @xs = ();
539                    my @ys = ();
540                    my @zs = ();
541                    my ($xv,$yv,$zv);
542                    while(@coords) {
543                        $xv = shift(@coords);
544                        $yv = shift(@coords);
545                        $zv = shift(@coords);
546                        $xv *= $magnification if $magX;
547                        $yv *= $magnification if $magY;
548                        $zv *= $magnification if $magZ;
549                        push @xs, ($xv*$T[3] + $yv*$T[6] + $zv*$T[9]  + $T[0]);
550                        push @ys, ($xv*$T[4] + $yv*$T[7] + $zv*$T[10] + $T[1]);
551                        push @zs, ($xv*$T[5] + $yv*$T[8] + $zv*$T[11] + $T[2]);
552                    }
553                    $elem->{'X'} = \@xs;
554                    $elem->{'Y'} = \@ys;
555                    $elem->{'Z'} = \@zs;
556                }
557            }
558        }
559        $i++;
560    }
561    if($xmin == $xmax) {
562        $xmin--;
563        $xmax++;
564    }
565    if($ymin == $ymax) {
566        $ymin--;
567        $ymax++;
568    }
569    if($zmin == $zmax) {
570        $zmin--;
571        $zmax++;
572    }
573    $xmin -= ($xmax - $xmin) / 6;
574    $xmax += ($xmax - $xmin) / 6;
575    $ymin -= ($xmax - $xmin) / 6;
576    $ymax += ($xmax - $xmin) / 6;
577    $zmin -= ($xmax - $xmin) / 6;
578    $zmax += ($xmax - $xmin) / 6;
579    %instr = ('x' => \@x, 'y' => \@y, z => \@z,
580              ori => \@ori, dis => \@dis, comp => \@comp,
581              xmin => $xmin, xmax => $xmax,
582              ymin => $ymin, ymax => $ymax,
583              zmin => $zmin, zmax => $zmax,
584              zoom_xmin => $xmin, zoom_xmax => $xmax,
585              zoom_ymin => $ymin, zoom_ymax => $ymax,
586              zoom_zmin => $zmin, zoom_zmax => $zmax,
587              zoom_tmin => 0, zoom_tmax => $tmax);
588    return %instr;
589}
590
591
592sub transform {
593    my ($comp, $x, $y, $z, $vx, $vy, $vz, $t, $ph1, $ph2) = @_;
594    if(!$comp) {
595        return ($x, $y, $z, $vx, $vy, $vz, $t, $ph1, $ph2);
596    } else {
597        my ($nx, $ny, $nz, $nvx, $nvy, $nvz, $nt, $nph1, $nph2);
598        my @T = @{$transformations{$comp}};
599        $x *= $magnification if $compdraw{$comp} && $compdraw{$comp}{'magX'};
600        $y *= $magnification if $compdraw{$comp} && $compdraw{$comp}{'magY'};
601        $z *= $magnification if $compdraw{$comp} && $compdraw{$comp}{'magZ'};
602        $nx = $x*$T[3] + $y*$T[6] + $z*$T[9] + $T[0];
603        $ny = $x*$T[4] + $y*$T[7] + $z*$T[10] + $T[1];
604        $nz = $x*$T[5] + $y*$T[8] + $z*$T[11] + $T[2];
605        $nvx = $vx*$T[3] + $vy*$T[6] + $vz*$T[9];
606        $nvy = $vx*$T[4] + $vy*$T[7] + $vz*$T[10];
607        $nvz = $vx*$T[5] + $vy*$T[8] + $vz*$T[11];
608        return ($nx, $ny, $nz, $nvx, $nvy, $nvz, $t, $ph1, $ph2);
609    }
610}
611
612
613sub get_inspect_pos {
614    my ($inspect, @comps) = @_;
615    return 0 unless $inspect;
616    my $i;
617    for($i = 0; $i < @comps; $i++) {
618        return $i if $comps[$i] eq $inspect;
619    }
620    die "Error: Inspected component $inspect not part of instrument $sim ?";
621}
622
623
624sub read_neutron {
625    my ($in) = @_;
626    my (@x, @y, @z, @vx, @vy, @vz, @t, @ph1, @ph2, @ncomp);
627    my ($st, $i);
628    my ($comp, $numcomp, $EndFlag);
629
630    $EndFlag = 0;
631    $st = 0;
632    $i = 0;
633    $numcomp = 0;
634    $EndFlag = 0;
635    my $dropit = 1;                # Flag to drop uninteresting neutron states.
636    while(<$in>) {
637        if($st == 0 && /^ENTER:/) {
638            # Neutron enters instrument.
639            $st = 1;
640        } elsif($st == 0 && /^STATE:/) {
641            # State after leaving - should probably be removed in McStas.
642            next;
643        } elsif($st == 1 && /^COMP:\s*"([a-zA-Z0-9_]+)"\s*$/) {
644            # Neutron enters component local coordinate system.
645            $comp = $1;
646            $numcomp++;
647            $dropit = 1;        # Drop the first state (entry point).
648        } elsif($st == 1 && (/^STATE:(.*)$/ || /^SCATTER:(.*)$/)) {
649            # Neutron state.
650            $dropit = 0, next if $dropit; # Skip entry point
651            ($x[$i], $y[$i], $z[$i],
652             $vx[$i], $vy[$i], $vz[$i],
653             $t[$i], $ph1[$i], $ph2[$i]) = split ",", $1;
654            ($x[$i], $y[$i], $z[$i],
655             $vx[$i], $vy[$i], $vz[$i],
656             $t[$i], $ph1[$i], $ph2[$i]) =
657                 transform($comp, $x[$i], $y[$i], $z[$i],
658                           $vx[$i], $vy[$i], $vz[$i],
659                           $t[$i], $ph1[$i], $ph2[$i]);
660            $ncomp[$i] = $comp;
661            if($TOF){
662                $t[$i] = 1000*$t[$i]; # Units of milli-seconds
663            }
664            $i++;
665        } elsif($st == 1 && /^ABSORB:/) {
666            # Neutron was absorbed.
667            next;                # No special action needed.
668        } elsif($st == 1 && /^LEAVE:/) {
669            # Neutron leaves instrument.
670            $st = 2;
671            last;
672        } elsif (/^Detector:/){
673          $st = 2;
674          if (! $MCSTAS::mcstas_config{'PLOTTER'} =~ /scriptfile/i) {
675            # Should only be done if finished, and not called with --save flag...
676            # Also, can only be done if tcl/tl available
677            if ($MCSTAS::mcstas_config{'TCLTK'} ne "no" && $EndFlag == 0) {
678              my $main = new MainWindow;
679              $main->Label(-text => "Simulation $sim ended."
680                          )->pack;
681              $main->Label(-text => 'Press Ok to terminate backend'
682                          )->pack;
683              $main->Button(-text => 'Ok',
684                            -command => sub{destroy $main}
685                           )->pack;
686              MainLoop;
687              $EndFlag = 1;
688            }
689          }
690          print;
691        } else {
692            # Pass on any output not meant for us.
693            print;
694        }
695    }
696    exit unless $st == 2;        # Stop when EOF seen before neutron data end.
697
698    my %neutron = ('x' => \@x, 'y' => \@y, z => \@z,
699                   vx => \@vx, vy => \@vy, vz => \@vz,
700                   t => \@t, ph1 => \@ph1, ph2 => \@ph2,
701                   comp => \@ncomp, numcomp => $numcomp, EndFlag => $EndFlag);
702    return %neutron
703}
704
705
706sub plot_components { # PGPLOT stuff only
707    my ($rx, $ry, $rori, $rdis, $axis1, $axis2) = @_;
708    my (@x, @y, @ori);
709    my ($i, $col);
710
711    @x = @$rx;
712    @y = @$ry;
713    @ori = @$rori;
714
715    PGPLOT::pgsci(2);
716    if ($TOF) {
717        my $zz;
718        my @ZZ;
719        my @TT;
720        $col = 4;
721        foreach $zz (@x) {
722            @ZZ = ($zz, $zz);
723            @TT = ($tmin, $tmax);
724            PGPLOT::pgsci($col++);
725            $col = 4 if $col > 15;
726            PGPLOT::pgline(2, \@TT, \@ZZ);
727        }
728    } else {
729        PGPLOT::pgline($#x + 1, \@x, \@y);
730        PGPLOT::pgpt($#x + 1, \@x, \@y, 20);
731        $col = 4;
732        for($i = 0; $i <= $#components; $i++) {
733            my $comp = $components[$i];
734            PGPLOT::pgsci($col++);
735            $col = 4 if $col > 15;
736            if($compdraw{$comp}) {
737                foreach $elem (@{$compdraw{$comp}{'elems'}}) {
738                    if($elem->{'type'} eq 'multiline') {
739                      PGPLOT::pgline($elem->{'count'}, $elem->{$axis1}, $elem->{$axis2});
740                    }
741                }
742            } else {
743                PGPLOT::pgsch(1.4);
744                  PGPLOT::pgpt(1, $x[$i], $y[$i], 26);
745              }
746        }
747    }
748}
749
750
751sub plot_neutron {
752    my ($x, $y, $z, $vx, $vy, $vz, $comp) = @_;
753    my ($i, $col, $oldcomp, $retval);
754    if ($MCSTAS::mcstas_config{'PLOTTER'} =~ /McStas|PGPLOT/i) {
755      # PGPLOT only
756      PGPLOT::pgsci(3);
757      PGPLOT::pgline(scalar(@$x), $x, $y);
758      # Show component entry/exit points in same colour as respective component.
759      # This should also be done w/ Matlab/Scilab...
760      $i = 0;
761      $col = 4;
762      while($i < scalar(@$x)) {
763        if(!defined($oldcomp) || $oldcomp cmp $comp->[$i]) {
764          $oldcomp = $comp->[$i];
765            PGPLOT::pgsci($col++);
766          $col = 4 if $col > 15;
767        }
768        # Exit point.
769        PGPLOT::pgpt(1, $x->[$i], $y->[$i], 17);
770        $i++;
771      }
772    } elsif ($MCSTAS::mcstas_config{'PLOTTER'} =~ /Matlab/i) {
773      # Matlab (split across multiple lines - otherwise sometimes
774      # crashes with component-rich instrs.)
775      # - Unless on Win32 where this will not work with the OLE connection
776      $retval=write_process("mcdisplay('Timeout');\n");
777      if ($Config{'osname'} eq 'MSWin32' && (!$file_output)) {
778        $retval=write_process("mcdisplay('PlotNeutron',[@$x],[@$y],[@$z]);\n");
779      } else {
780        $retval=write_process("mcdisplay('PlotNeutron',...\n");
781        $retval=write_process("[@$x],...\n");
782        $retval=write_process("[@$y],...\n");
783        $retval=write_process("[@$z]);\n");
784      }
785      return $retval;
786    } elsif ($MCSTAS::mcstas_config{'PLOTTER'} =~ /Scilab/i) {
787      # Scilab
788      $retval=write_process("x=[];\n");
789      $retval=write_process("y=[];\n");
790      $retval=write_process("z=[];\n");
791      $i=0;
792      while($i < scalar(@$x)) {
793        $tmp=$x->[$i];
794        $retval=write_process("x=[x $tmp];\n");
795        $tmp=$y->[$i];
796        $retval=write_process("y=[y $tmp];\n");
797        $tmp=$z->[$i];
798        $retval=write_process("z=[z $tmp];\n");
799        $i++;
800      }
801      $retval=write_process("PlotNeutron(x,y,z);\n");
802      return $retval;
803    } elsif ($MCSTAS::mcstas_config{'PLOTTER'} =~ /vrml/i) {
804      write_process("\nShape {\nappearance Appearance {\n\t".
805                    "material Material { emissiveColor 1 1 1\n transparency 0.7}}\n".
806                    "geometry IndexedLineSet {\ncoord Coordinate {\npoint [\n");
807      my $linelength="";
808      while($i < scalar(@$x)) {
809          my $xx,$yy,$zz;
810          # Needs swapping because of coordinate system in VRML instrument...
811          $xx = $y->[$i];
812          $yy = $z->[$i];
813          $zz = $x->[$i];
814          write_process("$xx $yy $zz, ");
815          $linelength="$linelength$i ";
816          $i++;
817      }
818      write_process("\n]}\ncoordIndex [$linelength -1,\n]}}");
819      return 1;
820    }
821}
822
823
824sub show_comp_names { # PGPLOT stuff only
825    my ($rinstr) = @_;
826    my %instr = %$rinstr;
827    my @comps = @{$instr{'comp'}};
828    my $count = @comps;
829    $count = 8 if $count < 8;
830    my $col = 4;
831    PGPLOT::pgsch(25/$count);
832    my $i;
833    for($i = 0; $i < @comps; $i++) {
834        PGPLOT::pgsci($col++);
835          $col = 4 if $col > 15;
836        if ($TOF) {
837            PGPLOT::pgmtxt('RV', 0.2, ($i+0.5)/$count, 0.0, $comps[$i]);
838        } else {
839            PGPLOT::pgmtxt('RV', 0.2, 1 - ($i+0.5)/$count, 0.0, $comps[$i]);
840        }
841    }
842}
843
844
845sub reset_zoom {
846    my ($rinstr, $vps) = @_;
847    $rinstr->{'zoom_xmin'} = $rinstr->{'xmin'};
848    $rinstr->{'zoom_xmax'} = $rinstr->{'xmax'};
849    $rinstr->{'zoom_ymin'} = $rinstr->{'ymin'};
850    $rinstr->{'zoom_ymax'} = $rinstr->{'ymax'};
851    $rinstr->{'zoom_zmin'} = $rinstr->{'zmin'};
852    $rinstr->{'zoom_zmax'} = $rinstr->{'zmax'};
853    $zooming = 0;
854}
855
856
857sub do_zoom {
858    my ($rinstr, $vps, $cx, $cy, $cx1, $cy1) = @_;
859    my ($tmp, $a1, $a2);
860    $tmp = $cx, $cx = $cx1, $cx1 = $tmp if $cx > $cx1;
861    $tmp = $cy, $cy = $cy1, $cy1 = $tmp if $cy > $cy1;
862
863    if($cx == $cx1 || $cy == $cy1) {
864        print STDERR "Warning: bad zoom area.\n";
865        return;
866    }
867    if($multi_view) {
868        # Convert from screen coordinates to world coordinates.
869        # First find which of the three views was choosen.
870        if($cx < 0 && $cy < 1) {
871            $cx = $cx + 1;
872            $cx1 = $cx1 + 1;
873            ($a1,$a2) = ("z", "y");
874        } elsif($cx < 0 && $cy >= 1) {
875            $cx = $cx + 1;
876            $cx1 = $cx1 + 1;
877            $cy = $cy - 1;
878            $cy1 = $cy1 - 1;
879            ($a1,$a2) = ("z", "x");
880        } elsif($cx >= 0 && $cy >= 1) {
881            $cy = $cy - 1;
882            $cy1 = $cy1 - 1;
883            ($a1,$a2) = ("x", "y");
884        } else {
885            print STDERR "Warning: bad zoom area.\n";
886            return;
887        }
888        my $idx = "$a1-$a2";
889        my $vpx0 = $vps->{$idx}{'VP'}[0];
890        my $vpdx = $vps->{$idx}{'VP'}[1] - $vpx0;
891        my $wx0 = $vps->{$idx}{'W'}[0];
892        my $wdx = $vps->{$idx}{'W'}[1] - $wx0;
893        my $vpy0 = $vps->{$idx}{'VP'}[2];
894        my $vpdy = $vps->{$idx}{'VP'}[3] - $vpy0;
895        my $wy0 = $vps->{$idx}{'W'}[2];
896        my $wdy = $vps->{$idx}{'W'}[3] - $wy0;
897        $cx = ($cx-$vpx0)/$vpdx*$wdx+$wx0;
898        $cx1 = ($cx1-$vpx0)/$vpdx*$wdx+$wx0;
899        $cy = ($cy-$vpy0)/$vpdy*$wdy+$wy0;
900        $cy1 = ($cy1-$vpy0)/$vpdy*$wdy+$wy0;
901    } else {
902        ($a1, $a2) = ("z","x");
903    }
904    $rinstr->{"zoom_${a1}min"} = $cx;
905    $rinstr->{"zoom_${a1}max"} = $cx1;
906    $rinstr->{"zoom_${a2}min"} = $cy;
907    $rinstr->{"zoom_${a2}max"} = $cy1;
908    $zooming = 1;
909    if ($TOF) {
910        $TOFINIT = 0;
911        $rinstr->{"zoom_tmin"} = $cx;
912        $rinstr->{"zoom_tmax"} = $cx1;
913        $rinstr->{"zoom_zmin"} = $cy;
914        $rinstr->{"zoom_zmax"} = $cy1;
915        $tmax=cx1;
916    }
917}
918
919sub write_process {
920  my ($command) = @_;
921  # $pid == 0 covers file output, Scilab/Matlab + Matlab on Win32 (OLE)
922  if (!$pid eq 0) {
923    (kill 0, $pid) || print STDERR "$plotter process terminated - ending...\n";
924    return 2;
925  }
926  if ($Config{'osname'} eq 'MSWin32' && $MCSTAS::mcstas_config{'PLOTTER'} =~ /Matlab/i && (!$file_output)) {
927    $ML->Execute($command);
928    my $err = Win32::OLE::LastError();
929    if ($err) {
930      print STDERR "Matlab terminated! - Exiting.\n";
931      return 2;
932    }
933  } else {
934    # Simply write data to pipe/file
935    print WRITER $command;
936    return 1;
937  }
938}
939
940sub vrml_setcolor       {
941    ($num,$max)=@_;# 0<=$num<$max
942        if($num % 2)
943    {$num=$num/2;}
944    else
945    {$num=$num/2+$max/2;}
946    my $num=sprintf("%d",$num);#floor
947        my $H=(6*$num)/$max;# 0<=Hue<6
948        my $iH=sprintf("%d",$H);#integer form 0 to 5
949        #Saturation and Value are fixed to 1
950        my $dH=$H-$iH;
951    my $R,$G,$B;
952    if($iH==0){$R=1    ;$G=$dH  ;$B=0    }
953    elsif($iH==1){$R=1-$dH;$G=1    ;$B=0    }
954    elsif($iH==2){$R=0    ;$G=1    ;$B=$dH  }
955    elsif($iH==3){$R=0    ;$G=1-$dH;$B=1    }
956    elsif($iH==4){$R=$dH  ;$G=0    ;$B=1    }
957    else{$R=1    ;$G=0    ;$B=1-$dH}
958    return "$R $G $B";
959}
960
961sub plot_instrument {
962    my ($noninteractive, $rinstr, $rneutron) = @_;
963    my %instr = %$rinstr;
964    my %neutron = %$rneutron;
965    my $retval;
966    # The following only relevant in PGPLOT mode
967    if ($MCSTAS::mcstas_config{'PLOTTER'} =~ /McStas|PGPLOT/i) {
968      my ($xmin, $xmax, $ymin, $ymax, $zmin, $zmax) =
969         ($instr{'zoom_xmin'}, $instr{'zoom_xmax'}, $instr{'zoom_ymin'},
970         $instr{'zoom_ymax'}, $instr{'zoom_zmin'}, $instr{'zoom_zmax'});
971      ($tmin, $tmax) = ($instr{'zoom_tmin'}, $instr{'zoom_tmax'});
972      my %vps;                        # Viewport/window setup.
973      my ($vpx1,$vpx2,$vpy1,$vpy2,$wx1,$wx2,$wy1,$wy2);
974
975      # PGPLOT::pgpage;        # start new page/panel
976      PGPLOT::pgbbuf;        # begin buffer batch output
977
978      # First show instrument from "above" (view in direction of y axis).
979
980      PGPLOT::pgsci(1);
981      PGPLOT::pgsch(1.4);
982      if ($TOF) {
983        if ($TOFINIT==0) {
984            PGPLOT::pgenv($tmin, $tmax, $zmin, $zmax, ($zooming ? 0 : 1), 0);
985            PGPLOT::pglab("TOF [ms]", "Z Axis [m]", "TOF diagram: $sim_cmd");
986            $TOFINIT=1;
987            show_comp_names($rinstr);
988        }
989        plot_components($instr{'z'}, $instr{'x'},$instr{'ori'}, $instr{'dis'},'T', 'Z');
990        plot_neutron($neutron{'t'}, $neutron{'z'}, $neutron{'y'},
991                     $neutron{'vz'}, $neutron{'vx'}, $neutron{'vy'},$neutron{'comp'});
992      } else {
993          if (!($keep) || ($keep && $PGINIT==0)) {
994              PGPLOT::pgenv($zmin, $zmax, $xmin, $xmax, ($zooming ? 0 : 1), 0);
995              PGPLOT::pglab("Z Axis [m]", "X Axis [m]", ($multi_view ? "Z-X view" : "Z-X view: $sim_cmd"));
996              show_comp_names($rinstr);
997              $PGINIT=1;
998          }
999        PGPLOT::pgsch(1.4);
1000        plot_components($instr{'z'}, $instr{'x'}, $instr{'ori'}, $instr{'dis'},
1001                      'Z', 'X');
1002        plot_neutron($neutron{'z'}, $neutron{'x'}, $neutron{'y'},
1003                     $neutron{'vz'}, $neutron{'vx'}, $neutron{'vy'},$neutron{'comp'});
1004      }
1005      if($multi_view) {
1006        # Remember viewport setup for Z-X view.
1007        PGPLOT::pgqvp(0, $vpx1, $vpx2, $vpy1, $vpy2);
1008        PGPLOT::pgqwin($wx1, $wx2, $wy1, $wy2);
1009        $vps{'z-x'} = {VP => [$vpx1,$vpx2,$vpy1,$vpy2],
1010                       W => [$wx1,$wx2,$wy1,$wy2]};
1011
1012        # Now show instrument viewed in direction of z axis.
1013        PGPLOT::pgsci(1);
1014        PGPLOT::pgsch(1.4);
1015        PGPLOT::pgenv($xmin, $xmax, $ymin, $ymax, ($zooming ? 0 : 1), 0);
1016        PGPLOT::pglab("X Axis [m]", "Y Axis [m]", "X-Y view");
1017        plot_components($instr{'x'}, $instr{'y'}, $instr{'ori'}, $instr{'dis'},
1018                        'X', 'Y');
1019        plot_neutron($neutron{'x'}, $neutron{'y'}, $neutron{'z'},
1020                     $neutron{'vx'}, $neutron{'vy'}, $neutron{'vz'} ,$neutron{'comp'});
1021        # Remember viewport setup for Z-X view.
1022        PGPLOT::pgqvp(0, $vpx1, $vpx2, $vpy1, $vpy2);
1023        PGPLOT::pgqwin($wx1, $wx2, $wy1, $wy2);
1024        $vps{'x-y'} = {VP => [$vpx1,$vpx2,$vpy1,$vpy2],
1025                       W => [$wx1,$wx2,$wy1,$wy2]};
1026
1027        # Now show instrument viewed in direction of x axis.
1028        PGPLOT::pgsci(1);
1029        PGPLOT::pgsch(1.4);
1030        PGPLOT::pgenv($zmin, $zmax, $ymin, $ymax, ($zooming ? 0 : 1), 0);
1031        PGPLOT::pglab("Z Axis [m]", "Y Axis [m]", "Z-Y view");
1032        plot_components($instr{'z'}, $instr{'y'}, $instr{'ori'}, $instr{'dis'},
1033                        'Z', 'Y');
1034        plot_neutron($neutron{'z'}, $neutron{'y'}, $neutron{'x'},
1035                     $neutron{'vz'}, $neutron{'vy'}, $neutron{'vx'}, $neutron{'comp'});
1036        # Remember viewport setup for Z-Y view.
1037        PGPLOT::pgqvp(0, $vpx1, $vpx2, $vpy1, $vpy2);
1038        PGPLOT::pgqwin($wx1, $wx2, $wy1, $wy2);
1039        $vps{'z-y'} = {VP => [$vpx1,$vpx2,$vpy1,$vpy2],
1040                       W => [$wx1,$wx2,$wy1,$wy2]};
1041
1042        # Set up viewport & window for mouse zoom.
1043        if ($multi_view) {
1044                PGPLOT::pgpanl(2,2);
1045                PGPLOT::pgsci(1);
1046                my $time=gmtime;
1047                PGPLOT::pgmtxt("t",0-1*1.2,0.0,0.0,"Date: $time");
1048                PGPLOT::pgmtxt("t",-2-1*1.2,0,0.0,"Simulation: ");
1049                PGPLOT::pgmtxt("t",-3-1*1.2,0.05,0.0,"$sim_cmd");
1050        }        # go to last panel
1051        else { PGPLOT::pgpanl(1,1); }                    # to activate full page
1052        PGPLOT::pgsvp(0,1,0,1);        # zoom for full page
1053        PGPLOT::pgswin(0,1,0,1);
1054      }
1055      PGPLOT::pgebuf;        # end buffer batch output
1056
1057      return 0 if ($noninteractive || $TOF || $keep);
1058
1059      # Now wait for a keypress in the graphics window.
1060      my ($cx, $cy, $cc);
1061      $cx = $cy = 0;
1062      PGPLOT::pgband(0, 0, 0, 0, $cx, $cy, $cc);
1063      if($cc =~ /[qQ]/) {
1064        return 2;                # Finished.
1065      } elsif($cc =~ /[pP]/) {        # B&W hardcopy.
1066        return 3;
1067      } elsif($cc =~ /[cC]/) {        # color hardcopy.
1068        return 4;
1069      } elsif($cc =~ /[gG]/) {        # GIF hardcopy.
1070        return 5;
1071      } elsif($cc =~ /[nN]/) {        # PNG hardcopy.
1072        return 6;
1073     } elsif($cc =~ /[hH]/) {        # Help
1074        print STDERR "McDisplay: q=quit, h=help\n";
1075        print STDERR "    output p=ps,   c=color ps, g=gif n=png\n";
1076        print STDERR "      zoom x=reset,z or d=zoom\n";
1077      } elsif($cc =~ /[zZdD]/) {        # Zoom.
1078        my ($cx1, $cy1, $cc1) = (0, 0, 0);
1079        PGPLOT::pgband(2,0,$cx,$cy,$cx1,$cy1,$cc1);
1080        do_zoom($rinstr, \%vps, $cx, $cy, $cx1, $cy1);
1081        return 1;
1082      } elsif($cc =~ /[xX]/) {        # Reset zoom.
1083              PGPLOT::pgpanl(1,1);
1084        PGPLOT::pgperas;
1085        if ($multi_view) {
1086                PGPLOT::pgpanl(1,2);
1087                PGPLOT::pgperas;
1088                PGPLOT::pgpanl(2,1);
1089                PGPLOT::pgperas;
1090                PGPLOT::pgpanl(2,2);
1091                PGPLOT::pgperas;
1092        }
1093        reset_zoom($rinstr, \%vps);
1094        return 1;
1095      }
1096    } else {
1097      # Leave further checks for plot_neutron
1098      $retval=plot_neutron($neutron{'z'}, $neutron{'x'}, $neutron{'y'},
1099                   $neutron{'vz'}, $neutron{'vx'}, $neutron{'vy'}, $neutron{'comp'});
1100      if ($retval==2) {
1101        return $retval;
1102      }
1103    }
1104    return 0;                        # Default: do not repeat this neutron.
1105}
1106
1107sub get_device { # PGPLOT stuff only
1108    my ($what) = @_;
1109    my $dev;
1110
1111    if (defined(&dev)) { $dev = dev($what); }
1112    else { $dev = PGPLOT::pgopen($what); }
1113    return $dev if $dev < 0;
1114    if($multi_view) {
1115    # We use a 2x2 display format to view the instrument from three angles.
1116        PGPLOT::pgsubp(2, 2);
1117    } else {
1118        # We use a 1x1 display for detail.
1119        PGPLOT::pgsubp(1, 1);
1120    }
1121    return $dev;
1122}
1123
1124
1125# Test the code.
1126
1127# Check command line arguments.
1128
1129undef $inspect;
1130undef $first;
1131undef $last;
1132undef $save;
1133undef $direct_output;
1134undef $sim_cmd;
1135undef $sim;
1136undef $TOF;
1137undef $tmax;
1138undef $keep;
1139undef $PGINIT;
1140undef $paramfile;
1141my $plotter;
1142undef $file_output;
1143my $int_mode=0; # interactive mode(0), non interactive (1)
1144my $i;
1145my $start_scilab=0;
1146my $show_help=0;
1147
1148$plotter = $MCSTAS::mcstas_config{'PLOTTER'};
1149
1150for($i = 0; $i < @ARGV; $i++) {
1151    if(($ARGV[$i] eq "-m") || ($ARGV[$i] eq "--multi")) {
1152        $multi_view = 1;
1153    } elsif($ARGV[$i] =~ /--help|-h$/) {
1154        $show_help=1;
1155    } elsif(($ARGV[$i] =~ /^-z([-0-9+.eE]+)$/) ||
1156            ($ARGV[$i] =~ /^--zoom=([-0-9+.eE]+)$/)) {
1157        $magnification = ($1 == 0 ? 1 : $1);
1158    } elsif(($ARGV[$i] eq "-gif") || ($ARGV[$i] eq "-ps") ||
1159            ($ARGV[$i] eq "-fig") || ($ARGV[$i] eq "-scg") ||
1160            ($ARGV[$i] eq "-psc") || ($ARGV[$i] eq "-png") || ($ARGV[$i] eq "-ppm")) {
1161        $direct_output = $ARGV[$i];
1162        $save = 1;
1163        $int_mode = 1;
1164    } elsif(($ARGV[$i] =~ /^-i([a-zA-Z0-9_]+)$/) ||
1165            ($ARGV[$i] =~ /^--inspect=([a-zA-Z0-9_]+)$/)) {
1166        $inspect = $1;
1167    } elsif($ARGV[$i] =~ /^--first=([a-zA-Z0-9_]+)$/) {
1168        $first = $1;
1169    } elsif($ARGV[$i] =~ /^--last=([a-zA-Z0-9_]+)$/) {
1170        $last = $1;
1171    } elsif($ARGV[$i] eq "--save") {
1172        $save = 1;
1173    } elsif(($ARGV[$i] =~ /^-p([a-zA-Z0-9_]+)$/) ||
1174              ($ARGV[$i] =~ /^--plotter=([a-zA-Z0-9_\"]+)$/) ||
1175              ($ARGV[$i] =~ /^--format=([a-zA-Z0-9_\"]+)$/)) {
1176        $plotter = $1;
1177   } elsif(($ARGV[$i] =~ /^-f([a-zA-Z0-9_\-\/\ \.\:\"]+)$/) ||
1178              ($ARGV[$i] =~ /^--file=([a-zA-Z0-9_\-\/\ \.\:]+)$/)) {
1179        $file_output = $1;
1180   } elsif($ARGV[$i] =~ /^--param=([a-zA-Z0-9_\ \"\.\-\:]+)$/) {
1181        $paramfile = $1;
1182   } elsif($ARGV[$i] eq "--TOF" || $ARGV[$i] eq "-T") {
1183       $TOF = 1;
1184   } elsif($ARGV[$i] eq "--keep" || $ARGV[$i] eq "-k") {
1185       $keep = 1;
1186       $int_mode = 0;
1187       $PGINIT = 0;
1188   } elsif($ARGV[$i] =~ /^--tmax=([0-9\.]+)$/) {
1189       $tmax=$1;
1190   } else {
1191        if (defined($sim_cmd)) { push @cmdline, $ARGV[$i]; }
1192        else {
1193          $sim_cmd = $ARGV[$i];
1194          $sim=$sim_cmd;
1195          # Remove trailing .out or .exe extension
1196          $sim=~ s|.out\Z||;
1197          $sim=~ s|.exe\Z||;
1198          $sim=~ s|.instr\Z||;
1199        }
1200   }
1201}
1202if ($show_help) { undef $sim_cmd; }
1203die "Usage: mcdisplay [-mzipfh][-gif|-ps|-psc] Instr.out [instr_options] params
1204 -h        --help            Show this help
1205 -m        --multi           Show the three instrument side views
1206 -T        --TOF             Special Time Of Flight acceptance diagram mode
1207           --tmax=TMAX       Maxiumum TOF [ms] (defaults to 50 ms)
1208 -zZF      --zoom=ZF         Show zoomed view by factor ZF
1209 -iCOMP    --inspect=COMP    Show only trajectories reaching component COMP
1210           --param=FILE     Read input parameters from parameter file
1211 -pPLOTTER --plotter=PLOTTER Output graphics using {PGPLOT,Scilab,Matlab}
1212 -fFNAME   --file=FNAME      Output graphics commands to file FNAME
1213                             (Only used when PLOTTER = {Scilab, Matlab})
1214           --first=COMP      First component to visualize {Scilab, Matlab}
1215           --last=COMP       Last component to visualize {Scilab, Matlab}
1216 -k        --keep            Plot all neutrons events together (Primarily for
1217                             use with PGPLOT and -psc / -gif etc.)
1218           --save            Output a Scilab/Matlab figure file and exit
1219                             (Filename is Instr.scf / Instr.fig). Figure
1220                             files are used by mcgui.pl for visualising the
1221                             instrument. With PGPLOT, --save is nonfunctional.
1222                             With VRML, --save disables spaw of VRML viewer.
1223 -gif|-ps|-psc               Export figure as gif/b&w ps/color ps and exit
1224 When using -ps -psc -gif, the program writes the hardcopy file and exits.
1225 SEE ALSO: mcstas, mcdoc, mcplot, mcrun, mcgui, mcresplot, mcstas2vitess
1226 DOC:      Please visit http://www.mcstas.org/\n"
1227 unless $sim_cmd;
1228
1229if($paramfile) {
1230    open(IN, "<$paramfile") || die "mcdisplay: Failed to open parameter file '$paramfile'";
1231    print "\nmcdisplay: Parameters specified using file \"$paramfile\"\n";
1232    while(<IN>) {
1233        my $p;
1234        for $p (split) {
1235            if($p =~ /^([A-Za-z0-9_]+)\=(.*)$/) {
1236                $parm = $1;
1237                $val = $2;
1238            } 
1239            @vals = split(',',$val);
1240            if (@vals>1) {
1241                $val = ($vals[0]+$vals[1])/2;
1242                print "Parameter $parm: Substituting interval $vals[0],$vals[1] to $val\n";
1243                push @cmdline, "$parm=$val";
1244            } else {
1245                push @cmdline, $p; 
1246            }
1247        }
1248    }
1249    print "\n";
1250}
1251
1252if ($sim_cmd =~ m'\.instr$') # recompile .instr if needed
1253{ my @ccopts=();
1254  ($sim_cmd, undef) = get_out_file($sim_cmd, 0, @ccopts); }
1255
1256
1257# Check value of $plotter and $file_output variables, set
1258# $MCSTAS::mcstas_config{'PLOTTER'} with scriptfile keyword, always set for VRML
1259if ($file_output || $plotter =~ /VRML/i) { $plotter .= "_scriptfile"; }
1260
1261if ($TOF) {
1262    $TOFINIT=0;
1263    if(! ($tmax)) {
1264        $tmax=50;
1265    }
1266    if (!($plotter =~ /McStas|PGPLOT/i)) {
1267        print STDERR "\n***************************************\n";
1268        print STDERR "TOF only possible using plotter PGPLOT\nSelecting PGPLOT";
1269        print STDERR "\n***************************************\n";
1270        $plotter="PGPLOT";
1271    }
1272}
1273
1274if ($plotter =~ /Scilab/i && not $plotter =~ /scriptfile/i) {
1275  # Check for Win32, only file_output possible :(
1276  if ($Config{'osname'} eq MSWin32) {
1277    $file_output="mcdisplay_commands.sci";
1278    print STDERR "\n******************************************************\n";
1279    print STDERR "Sorry, Scilab only possible using file output on Win32\n\n";
1280    print STDERR "Outputting to file $file_output\n";
1281    print STDERR "******************************************************\n\n";
1282    print STDERR "When done, I will start a scilab for you to view the file...\n";
1283    $start_scilab=1;
1284    $plotter="Scilab_scriptfile";
1285  }
1286}
1287
1288if ($plotter =~ /scriptfile/i && not $file_output) {
1289  $file_output="mcdisplay_commands";
1290  if ($plotter =~ /Matlab/i) { $file_output .=".m"; }
1291  elsif ($plotter =~ /Scilab/i) { $file_output .=".sci"; }
1292  elsif ($plotter =~ /VRML/i) { $file_output .=".wrl"; }
1293  print STDERR "Outputting to file $file_output\n";
1294}
1295
1296# Final PLOTTER check, is PGPLOT wanted but not possible?
1297# - Ask user to rerun / set other default
1298if ($plotter =~ /McStas|PGPLOT/i && $MCSTAS::mcstas_config{'PGPLOT'} eq "no") {
1299  print STDERR "\n******************************************************\n";
1300  print STDERR "Default / selected PLOTTER is PGPLOT - Problems:\n\n";
1301  print STDERR "PGPLOT.pm not found on Perl \@INC path\n\nSolution:\n\n";
1302  print STDERR "1) Install pgplot + pgperl packages (Unix/Linux/Cygwin) \n";
1303  print STDERR "2) Rerun mcdisplay with -p/--plotter set to Scilab/Matlab \n";
1304  print STDERR "3) Modify $MCSTAS::perl_dir/mcstas_config.perl\n";
1305  print STDERR "   to set a different default plotter\n\n";
1306  print STDERR "******************************************************\n\n";
1307  die "mcdisplay: PGPLOT problems...\n";
1308}
1309
1310$MCSTAS::mcstas_config{'PLOTTER'} = $plotter;
1311
1312my $pg_devname = $MCSTAS::mcstas_config{'PGDEV'};
1313# Only set up PGPLOT stuff if needed
1314if ($plotter =~ /McStas|PGPLOT/i) { # PGPLOT is plotter!
1315  if ($int_mode == 1)
1316    {
1317      my $ext  = "ps";
1318      my $type = "ps";
1319      if($direct_output eq "-gif") { $ext="gif"; $type="gif"; }
1320      elsif($direct_output eq "-png") { $ext="png"; $type="png"; }
1321      elsif($direct_output eq "-psc") { $type="cps"; }
1322      $pg_devname = "$sim_cmd.$ext/$type";
1323    }
1324  $global_device = get_device($pg_devname);
1325  if($global_device < 0) {
1326    print STDERR "mcdisplay: Failed to open PGPLOT device $pg_devname\n";
1327    exit 1;
1328  }
1329  my $seq = "";                        # Sequence number for multiple hardcopy
1330  PGPLOT::pgask(0);
1331} elsif ($plotter =~ /vrml/i) {
1332  # VRML always with file handle. VRML browser spawned after write of file.
1333  print "Opening file ...\n";
1334  open(WRITER, "> $file_output");
1335  $pid=0;
1336  # Other VRML init stuff:
1337  my $nbcomp1=0,$nbcomp2=0;
1338  my %transforms;
1339  my %compheaders;
1340} elsif ($plotter =~ /Matlab/i && $plotter =~ /scriptfile/i) {
1341  # Matlab w/FILE is plotter - open a file handle
1342  open(WRITER, "> $file_output");
1343  $pid=0;
1344} elsif ($plotter =~ /Matlab/i) {
1345  # Matlab is plotter - open a pipe / OLE connection
1346  if ($Config{'osname'} eq 'MSWin32') {
1347    $pid=0;
1348    $ML = Win32::OLE->new('Matlab.Application') || die "mcdisplay: Could not start Matlab\n";
1349  } else {
1350    my $cmd = "$MCSTAS::mcstas_config{'MATLAB'} $MCSTAS::mcstas_config{'MATLAB_COMMAND'} > /dev/null";
1351    $pid=open2(READER,WRITER, $cmd) || die "mcdisplay: Could not start Matlab\n$cmd\n";
1352  }
1353  print STDERR "Opened up pipe to matlab - pid is $pid\n";
1354  print STDERR "Building Matlab INSTRUMENT struct, please have patience...\n";
1355} elsif ($plotter =~ /Scilab/i && $plotter =~ /scriptfile/i) {
1356  # Scilab w/FILE is plotter - open a file handle
1357  open(WRITER, "> $file_output");
1358  $pid=0;
1359} elsif ($plotter =~ /Scilab/i) {
1360  # Scilab is plotter - open a pipe
1361  my $cmd = "$MCSTAS::mcstas_config{'SCILAB'} $MCSTAS::mcstas_config{'SCILAB_COMMAND'} > /dev/null";
1362  $pid=open2(READER,WRITER, $cmd) || die "Could not start Scilab\n$cmd\n";
1363  print STDERR "Opened up pipe to scilab - pid is $pid\n";
1364  print STDERR "Building Scilab INSTRUMENT struct, please have patience...\n";
1365}
1366
1367my ($numcomp, %neutron, %instr);
1368
1369$args = join(" ", @cmdline);
1370$cmdline = "$sim_cmd --trace --no-output-files $args";
1371printf STDERR "Starting simulation '$cmdline' ...\n";
1372open(IN, "$cmdline |") || die "mcdisplay: Could not run simulation\n";
1373
1374$numcomp = read_instrument(IN);
1375$inspect_pos = get_inspect_pos($inspect, @components);
1376%instr = make_instrument;
1377
1378if ($int_mode == 0 && $plotter =~ /McStas|PGPLOT/i)  { printf STDERR "McDisplay/PGPLOT: Press H key for help.\n"; }
1379if ($plotter =~ /Matlab/i && not $plotter =~ /scriptfile/i) { print STDERR "Matlab INSTRUMENT done, starting gui....\n"; }
1380if ($plotter =~ /Scilab/i && not $plotter =~ /scriptfile/i) { print STDERR "Scilab INSTRUMENT done, starting gui....\n"; }
1381
1382while(!eof(IN)) {
1383    %neutron = read_neutron(IN);
1384    if ($start_scilab == 1) {
1385    # This only happens on Win32 (runscilab.exe), and we know the filename too...
1386        my $runscilab = "$MCSTAS::mcstas_config{'SCILAB'}";
1387        my $pid = open(SCILAB,"$runscilab -f mcdisplay_commands.sci|");
1388        while(!eof(SCILAB)) {
1389            # Do nothing...
1390        }
1391        $start_scilab = 0;
1392    }
1393    next if $neutron{'numcomp'} <= $inspect_pos;
1394
1395    my $ret;
1396    do {
1397        $ret = plot_instrument($int_mode, \%instr, \%neutron);
1398        if ($plotter =~ /McStas|PGPLOT/i) {
1399          if (! ($TOF) && ! ($keep)) {
1400                if ($int_mode == 1) { $ret =2; print STDERR "Wrote \"$pg_devname\"\n"; }
1401          }
1402          if($ret == 3 || $ret == 4 || $ret == 5 || $ret == 6) {
1403            my $ext="ps";
1404            my $type = $ret == 3 ? "ps" : "cps";
1405            if($ret == 5) { $type = "gif"; $ext="gif"; }
1406            if($ret == 6) { $type = "png"; $ext="png"; }
1407            my $tmp_pg_devname = "$sim_cmd$seq.$ext/$type";
1408            my $tmpdev=0;
1409            $tmpdev = get_device($tmp_pg_devname);
1410            if($tmpdev < 0) {
1411              print STDERR "mcdisplay: Warning: could not open PGPLOT output \"$tmp_pg_devname\" for hardcopy output\n";
1412            } else {
1413              plot_instrument(1, \%instr, \%neutron);
1414              print STDERR "Wrote \"$tmp_pg_devname\"\n";
1415              ++$seq;
1416            }
1417            if (defined(&close_window)) {
1418                close_window($tmpdev);
1419            } else {
1420                PGPLOT::pgclos();
1421            }
1422            PGPLOT::pgslct($global_device);
1423#           PGPLOT::pgask(0);
1424          }
1425        }
1426    } while($ret != 0 && $ret != 2);
1427    last if $ret == 2;
1428  }
1429close(IN);
1430if ($plotter =~ /VRML/i && !($MCSTAS::mcstas_config{'VRMLVIEW'} eq "no")) {
1431    if (-e $file_output) {
1432        if (! $save) {
1433            print STDERR "Spawning $MCSTAS::mcstas_config{'VRMLVIEW'} to view $file_output\n";
1434            open(VRML, "$MCSTAS::mcstas_config{'VRMLVIEW'} $file_output|");
1435        }
1436    }
1437}
1438
1439# Properly close any open files etc.
1440if ($plotter =~ /McStas|PGPLOT/i) {
1441  if (defined(&close_window)) { close_window(); }
1442  else {  PGPLOT::pgclos(); }
1443} else {
1444  close(WRITER);
1445}
1446
Note: See TracBrowser for help on using the browser.