#!/usr/bin/perl
#
# ./shower2pv airesoutput outfile steptime boxheight numsteps 
#
# ./shower2pv airesoutput [-o outfile] [-t steptime] [-b boxheight] [-n numsteps] [-static] [-lines/-thicklines] [-ground groundpicturename] [-groundsize lengthInMetersOfSideOfGroundSquare]
#
# 
# showerfile is as output by AiresQ (using "./AiresQ < example.inp > showerfile")
# This creates outfile.cf, and other Partiview-format files of the form outfile*
# 
# steptime is the number of nanoseconds that are considered to be a single time unit 
# boxheight is the height, in meters, of the box (0 if you dont want a box) that is 5km x 5km on the ground
# numsteps is the number of time units you want (useful if you're not interested in particles after some point)

# Defaults: only the first argument (airesoutput) has to be specified.
#   outfile : assumed to equal airesoutput if not specified.
#   steptime:  200ns if not specified
#   boxheight: set to 0 (i.e. no box) if not specified
#   numsteps : set to all if not specified
#   ground: ground.sgi
#   groundsize: 8000 meters
#   using dots

# This file works by creating a series of files tmp_d.txt  for all integers d from 0
#  to the maximum number of timesteps. Then these files are combined by calling the *nix
#  "cat" command repeatedly. For some reason I only wanted to cat two files at a time,
#  so it creates two files TMP0 and TMP1 that are cumulative collections of the tmp_d.txt files.
#
# The temporary files are deleted at the end
#
# Dinoj Surendran (dinoj@cs.uchicago.edu), Dec 2004 (earlier versions made in Aug 2004)
# See http://astro.uchicago.edu/cosmus/projects/aires
# Modified 5 Jan 2005, made all cmap outputs same, fixed bug with allocating colors for missing particles 
#   ... apparently aires generates undocumented particles (eg. particle 142 during iron showers)
#
# Corrected by Mark SubbaRao of COSMUS in Nov/December 2005
# 
# You need to name a file called ground.sgi

use bytes;
use POSIX qw(ceil floor);

my $airesfile = $ARGV[0];
my $outfile = $ARGV[0];
my $steptime = 200;
my $boxheight = 0;
my $numsteps = 0;
my $USELINES = 0;          # equals 1 for thin lines, 2 for thick lines
my $LASTFRAMEONLY = 0;
my $VERBOSE = 0;
my $i = 1;
my $ground = "ground";
my $groundsize = 8000;   # length of side, in meters, of square representing the ground (so the default assumption here is that it's 8km.)
my $toffset = 0;  # timestep at which the shower starts

while ($i <= $#ARGV)
{
#    print "argv[$i]=[", $ARGV[$i], "]\n";

    if ($ARGV[$i] eq "-t")
    {
	$i++;
	if ($i <= $#ARGV)
	{$steptime = $ARGV[$i];}
    }
    elsif ($ARGV[$i] eq "-o")
    {
	$i++;
	if ($i <= $#ARGV)
	{$outfile = $ARGV[$i];}
    }
    elsif ($ARGV[$i] eq "-b")
    {
	$i++;
	if ($i <= $#ARGV)
	{ $boxheight = $ARGV[$i];}
    }
    elsif ($ARGV[$i] eq "-n")
    {
	$i++;
	if ($i <= $#ARGV)
	{$numsteps = $ARGV[$i];}
    }
    elsif ($ARGV[$i] eq "-to")
    {
	$i++;
	if ($i <= $#ARGV)
	{$toffset = $ARGV[$i];}
    }
    elsif ($ARGV[$i] eq "-lines")
    {
	$USELINES = 1;
    }
    elsif ($ARGV[$i] eq "-thicklines")
    {
	$USELINES = 2;
    }
    elsif ($ARGV[$i] eq "-current")
    {
	$USELINES = 3;
    }
    elsif ($ARGV[$i] eq "-static")
    {
	$LASTFRAMEONLY = 1;
    }
    elsif ($ARGV[$i] eq "-verbose")
    {
	$VERBOSE = 1;
    }
    elsif ($ARGV[$i] eq "-ground")
    {
	$i++;
	if ($i <= $#ARGV)
	{$ground = $ARGV[$i];}
    }
    elsif ($ARGV[$i] eq "-groundsize")
    {
	$i++;
	if ($i <= $#ARGV)
	{$groundsize = $ARGV[$i];}
    }
    $i++;
}
if ($VERBOSE)
{
    print "steptime = $steptime\n";
    print "outfile = $outfile\n";
    print "boxheight = $boxheight\n";
    print "numsteps = $numsteps\n";
    print "uselines = $USELINES\n";
    print "static = $LASTFRAMEONLY\n";
}

#print "[ground = $ground]\n";
if ($ground !~ /\./)
  {$ground = $ground.".sgi";}
#print "[ground = $ground]\n";

$NUMPARTICLESIZES = 10;
$USELOGENERGY = 1;
$NUMREPORT = 5000; # report once every <this many> particles are processed

$BUTTONTHRESH = 1; #Create buttons only for points with at least BUTTONTHRESH appearances # not implemented yet

open (IN,$airesfile);
$showerfile = "shower_".$outfile;
open (SHOWER,">".$showerfile);
print SHOWER "% particletype energy timestart timeend x_start y_start z_start x_end y_end z_end\n";

my %typestats;    # number of particles of each type
my %typenames;

 # add colorname e.g. $colname{"1"} = "green" so cmap file is more interpretable

$typenames{"1"}="gamma";      $col{"1"} = "1.0 0.5 0";      # gamma is orange    
#$typenames{"1"}="gamma";      $col{"1"} = "0 1 0";      # gamma is green    

#$typenames{"2"}="positron";   $col{"2"} = "1 0 1";        # positron (e+) is purple
$typenames{"2"}="positron";   $col{"2"} = "0.5 0 0.7";        # positron (e+) is purple

#$typenames{"-2"}="electron";  $col{"-2"} = "0 1 1";        # electron (e-) is cyan
$typenames{"-2"}="electron";  $col{"-2"} = "0.7 0 0.5";        # electron (e-) is magenta

$typenames{"3"}="muonP";      $col{"3"} = "0.9 0.1 0.0 ";       # muon+ is red
$typenames{"-3"}="muonN";     $col{"-3"}= "1 0 0 ";       # muon- is red
$typenames{"4"}="tauP";       $col{"4"} = "1 1 1 ";       # tau+ is white
$typenames{"-4"}="tauN";      $col{"-4"}= "1 1 1 ";       # tau- is white
$typenames{"6"}="nuE";        $col{"6"} = "1 1 1 ";       # nu(e) is white
$typenames{"-6"}="nubarE";    $col{"-6"}= "1 1 1 ";       # nubar(e) is white
$typenames{"7"}="nuM";        $col{"7"}= "1 1 1 ";        # nu(m) is white
$typenames{"-7"}="nubarM";    $col{"-7"}= "1 1 1 ";       # nubar(m) is white
$typenames{"8"}="nuT";        $col{"8"}= "1 1 1 ";        # nu(t) is white
$typenames{"-8"}="nubarT";    $col{"-8"}= "1 1 1 ";       # nubar(t) is white
$typenames{"10"}="pi0";       $col{"10"}= "0 1 0 ";       # pi0 is green 
$typenames{"11"}="piP";       $col{"11"}= "0 1 0 ";       # pi+ is green 
$typenames{"-11"}="piN";      $col{"-11"}= "0 1 0 ";      # pi- is green 
$typenames{"12"}="K0S";       $col{"12"}= "1 1 1 ";       # K0S is white 
$typenames{"13"}="K0L";       $col{"13"}= "1 1 1 ";       # K0L is white 
$typenames{"14"}="KP";        $col{"14"}= "1 1 1 ";       # K+ is white 
$typenames{"-14"}="KN";       $col{"-14"}= "1 1 1 ";      # K- is white 
$typenames{"15"}="eta";       $col{"15"}= "1 1 1 ";       # eta is white 
$typenames{"20"}="lambdaA";   $col{"20"}= "1 1 1 ";       # lambdaA is white 
$typenames{"-20"}="lambdaB";  $col{"-20"}= "1 1 1 ";      # lambdaB is white 
$typenames{"30"}="neutron";   $col{"30"}= "1 1 0 ";       # neutron is yellow
$typenames{"-30"}="antineutron";$col{"-30"}= "1 0.8 0 ";    # antineutron is darker yellow
$typenames{"31"}="proton";    $col{"31"}= "0 0.0 1.0 ";       # proton is blue
$typenames{"-31"}="antiproton";$col{"-31"}= "0 0.1 1.0 ";      # antiproton is teal
$typenames{"0"}="other";   $col{"0"} = "0.5 0.5 0.5"; # any other particle is grey

foreach $t (keys (%typenames))
{
    $typestats{$t} = 0;
}

my $minenergy = 1e50;
my $maxenergy =-1e50;

my $mintime = 1e50;
my $maxtime =-1e50;

my $minlifetime = 1e50;
my $maxlifetime =-1e50;

my @energys;

# lines are of four types. 
# 1) PSTART : info about the start of a particle
# 2) PEND : info about the end of a particle
# 3) ------- i.e. dividing line
# 4) other lines - ignore them
 
$seenfirstofpair = 0;

while ($line = <IN>)
{
    chomp $line;
    $PSTART = 0; $PEND = 0; $DIVIDER = 0;
    @tmp = split /[\t\n\r ]+/, $line;

    if (($#tmp == 1) & ($tmp[1] =~ /---/))
    {
	$DIVIDER = 1;
    }
    elsif ($#tmp >= 6)
    {
	$type = $tmp[1];
	$x = $tmp[2];
	$y = $tmp[3];
	$z = $tmp[4];
	$engy = $tmp[5];
	$time = $tmp[6];

	if (isnumeric ($type) & (isnumeric ($x)))
	{
	    if (($#tmp >= 7) & ($tmp[7] =~ "Start"))
	    {
		$PSTART = 1;
	    }
	    else
	    {
		$PEND = 1;
	    }
	}
    }

    if (($PEND == 0) & ($PSTART == 1))
    {
	$sx = $x; $sy = $y; $sz = $z;
	$senergy = $engy; $stime = $time;
	$seenfirstofpair = 1;
    }
    elsif (($PEND == 1) & ($seenfirstofpair == 1))
    {
	$seenfirstofpair = 0;    ################# bug removed?

	$ex = $x; $ey = $y; $ez = $z;
	$eenergy = $energy; $etime = $time;
	$energy = 0.5*($senergy+$eenergy);
	$lifetime = $etime-$stime;

	if ($energy > $maxenergy) {$maxenergy = $energy;}
	if ($energy < $minenergy) {$minenergy = $energy;}
	if ($stime > $maxtime) {$maxtime = $stime;}
	if ($stime < $mintime) {$mintime = $stime;}
	if ($etime > $maxtime) {$maxtime = $etime;}
	if ($etime < $mintime) {$mintime = $etime;}
	if ($lifetime > $maxlifetime) {$maxlifetime = $lifetime;}
	if ($lifetime < $minlifetime) {$minlifetime = $lifetime;}

	if (($sx != 0) & ($senergy != 0))
	{
	    #print out info about this particle
	    printf SHOWER "%d %0.10f %0.3f %0.3f %0.5f %0.5f %0.5f %0.5f %0.5f %0.5f \n",$type, $energy, $stime, $etime, $sx, $sy, $sz, $ex, $ey, $ez;

	    $typestats{$type} = $typestats{$type} + 1;
	}
    }
    else
    {
	$sx = 0; $sy = 0; $sz = 0; 
	$ex = 0; $ey = 0; $ez = 0; 
	$senergy = 0; $stime = 0;
	$eenergy = 0; $etime = 0;
	$seenfirstofpair = 0;
    }
}

# dump type stats # not really required
close (SHOWER);
foreach $t (keys (%typestats))
{
    if (! exists ($typenames{$t}))
    {
	$typestats{0} = $typestats{0} + $typestats{$t};
    }
}

open (SUMMARY, ">summary_".$outfile);

printf SUMMARY "Energy values range from %0.5f to %0.5f\n", $minenergy, $maxenergy;
printf SUMMARY "Time values range from %0.5f to %0.5f\n", $mintime, $maxtime;
printf SUMMARY "Particle lifetime values range from %0.5f to %0.5f\n", $minlifetime, $maxlifetime;

@blah = keys (%typestats);

$cmapfile = $outfile.".cmap";
open (CMAP,">".$cmapfile);
printf CMAP "%d\n", $#blah+2;
print CMAP "1 1 1\n";
$colindex = {};            # first colindex should be 1
$count = 0;
foreach $t (keys (%typenames))
{
    printf SUMMARY "$t: %s : %d\n", $typenames{$t}, $typestats{$t};
    print CMAP $col{$t}, "# ", $typenames{$t}, " (type=", $t, ", number=", $typestats{$t},")\n";
    $count++;
    $colindex{$t} = $count;
}

close (SUMMARY);

######### Create buttons only for points with at least BUTTONTHRESH appearances ##### not implemented yet ####

if (0 >= $minenergy) {$minenergy = 1e-10;}

$minlogenergy = log($minenergy);

open (SHOWER, $showerfile) or die ("$showerfile not found\n");

$maxtimebinned=POSIX::ceil (($maxtime-$mintime)/$steptime);

$maxtimebinned;

printf ("\nIn this shower, the earliest time is at %d nanoseconds and the latest at %d ns.\n", $mintime, $maxtime);
printf ("  You want each time step to be %0.2f ns, so there will be %d time steps in this shower.\n", $steptime, $maxtimebinned+1);
if (($numsteps > 0) & ($numsteps < $maxtimebinned))
{
    $maxtimebinned = $numsteps;
    printf ("  Of these, only the first %d time steps will be used.\n", $maxtimebinned);
}
printf ("Type 'step %d' when in partiview to see the whole shower)\n\n", $maxtimebinned);

my @tmpfiles;

if ($LASTFRAMEONLY == 0)
{
    for ($i=0; $i<=$maxtimebinned; $i++)
    {
        $j=$i+$toffest;
	local *FILE;
	open(FILE, ">tmp_$i.txt") || die;
	push(@tmpfiles, *FILE);
	
	open(HED, ">header_$i.txt") || die;
	print HED "\ndatatime $j\n\n";
	close(HED);
    }

    open(HED, ">header.txt") ;
    print HED "datavar 0 colindex\ndatavar 1 energy\ndatavar 2 type\ntexturevar 3\n\n";
    close(HED);
}
elsif ($LASTFRAMEONLY == 1)
{
    open(SPECK, ">".$outfile.".speck") || die;
    print SPECK "datavar 0 colindex\ndatavar 1 energy\ndatavar 2 type\ntexturevar 3\n\n";
}

# print "[$maxtimebinned]  [$minenergy] [$maxenergy] [$mintime] [$maxtime] [$minlifetime] [$maxlifetime]\n";

printf "Number of particles processed:";
$count = 0;
while ($line = <SHOWER>)
{
    next if ($line =~ "%");
    ($type, $energy, $stime, $etime, $sx, $sy, $sz, $ex, $ey, $ez) = split /[\t\n\r ]+/, $line;

    $color = 0; 
    if (! exists ($colindex{$type}))
    {
	if (1 == $typestats{$type})
	{
	    print "\n  weird... a particle of unknown type ($type) occurred once.\n";
	}
	else
	{
	    print "\n  weird... a particle of unknown type ($type) occurred ",$typestats{$type}," times. \n";
	}
	$type = 0;
    }
    $color = $colindex{$type};


    if (1 == $USELOGENERGY)
    {
	$size = $NUMPARTICLESIZES * (log($energy)-$minlogenergy) / (log($maxenergy)-$minlogenergy);
    }
    else
    {
    	$size = $NUMPARTICLESIZES * (($energy)-($minenergy)) / (($maxenergy)-($minenergy));
    }
    $size=1 if (0==$size);

    # This particle appears between timesteps $st and $et inclusive

    $st = POSIX::ceil (($stime-$mintime) / $steptime);
    $et = POSIX::ceil (($etime-$mintime) / $steptime);

    if (($LASTFRAMEONLY == 1) & ($USELINES >= 1))
    {
	if ($USELINES == 1)
	{
		printf SPECK thinline ($sx,$sy,$sz,$ex,$ey,$ez,$size,$color), "\n";
	}
	elsif ($USELINES == 1)
	{
		printf SPECK thickline ($sx,$sy,$sz,$ex,$ey,$ez,$size,$color), "\n";
	}
    }
    else
    {
	$ft = min($maxtimebinned,$et);
	for ($t=$st; $t<=$ft; $t++)
	{
	    # Write it in $et-$st+1 files, assuming it moves with constant speed (that of light) during its life
	    $fraclife = 1.0*($t-$st+1)/($et-$st+1);   # fraction of life it's in at this time $t
	    
	    $x = $sx + $fraclife*($ex-$sx);
	    $y = $sy + $fraclife*($ey-$sy);
	    $z = $sz + $fraclife*($ez-$sz);
	    
	    if ($LASTFRAMEONLY == 0)
	    {
		$FILE = $tmpfiles[$t];
		if ($USELINES == 0)
		{
		    printf $FILE "%0.1f %0.1f %0.1f %d %d %d 1 #\n", $x,$y,$z,$color,$size,$type;
		}
		elsif ($USELINES >= 1)
		{
	            $s = "";
		    if ($USELINES == 1)
		    {
		       $s = thinline ($sx,$sy,$sz,$x,$y,$z,$size,$color);
	            }
	            elsif ($USELINES == 2)
		    {
		       $s = thickline ($sx,$sy,$sz,$x,$y,$z,$size,$color);
		    }	
		    printf $FILE "$s\n";
		}
	    }
	    elsif ($LASTFRAMEONLY == 1)  # must have USELINES == 0
	    {
		print SPECK "%0.1f %0.1f %0.1f %d %d %d 1 #\n", $x,$y,$z,$color,$size,$type;
	    }
	}
	if ($USELINES >= 1)
	{
	    $s = "";
	    if ($USELINES == 1)
	    {
	      $s = thinline ($sx,$sy,$sz,$ex,$ey,$ez,$size,$color);
	    }
	    elsif ($USELINES == 2)
	    {
	      $s = thickline ($sx,$sy,$sz,$ex,$ey,$ez,$size,$color);
	    }
	    for ($t=$ft+1; $t <= $maxtimebinned; $t++)
	    {
		$FILE = $tmpfiles[$t];
		printf $FILE "$s\n";
	    }
	}
    }

    printf ("...%d particles processed\n", $count) if (0==(++$count % $NUMREPORT)) 
}

printf "...%d particles processed\n\n", $count;

close (SHOWER);


if ($LASTFRAMEONLY == 1)
{
    close(SPECK);
}
else
{
#    if ($USELINES == 1)
#    {
#	$cmd = "cat header.txt ";
#	for ($i=0; $i<$maxtimebinned; $t++)
#	{
#	    $cmd = $cmd." header_$t.txt tmp_$t.txt ";
#	}
#	$cmd = $cmd. " > ".$outfile.".speck";
#
#	print "$cmd\n";
#
#	system ("$cmd");
#	system ("rm tmp*.txt header*.txt");
#    }
#    else
#    {
	open (TMP,">TMP0");
	print TMP " ";
	close (TMP);
	
	open (TMP,">TMP1");
	print TMP " ";
	close (TMP);
	for ($t=0; $t<=$maxtimebinned; $t++)
	{
	    if (0 == mod($t,2))
	    {
		$cmd = "cat TMP0 header_$t.txt ";
	    }
	    else
	    {
		$cmd = "cat TMP1 header_$t.txt ";
	    }
	

	    if ($USELINES == 0)
	    {
		for ($s=0; $s <= $t; $s++)
  		  {$cmd = $cmd." tmp_$s.txt ";} 
	    }
	    else
	    {
  		$cmd = $cmd." tmp_$t.txt ";
	    }

	    if (0 == mod($t,2))
	    {
		$cmd = $cmd . "> TMP1";
	    }
	    else
	    {
		$cmd = $cmd . "> TMP0";
	    }
	    
	    system ("$cmd");
	    printf ("%d of %d timesteps processed \n",$t, $maxtimebinned);
	}
    
	system ("rm tmp*.txt header_*.txt");
    
	if (0 == mod($maxtimebinned,2))
	{
	    $cmd = "cat header.txt TMP1 > ".$outfile.".speck";
	}
	else
	{
	    $cmd = "cat header.txt TMP0 > ".$outfile.".speck";
	}
	system ("$cmd");

	system ("rm TMP* he*.txt");
#    }
}

open (CF,">".$outfile.".cf");

printf CF "filepath +:.:./data\n\nobject g1=shower\n  include %s.speck\neval cmap %s.cmap\neval color colindex\neval lum energy\neval points on\neval poly off\neval ptsize 1 10\neval psize 1e8\neval slum 100\n\n", $outfile, $outfile;

printf CF "object g2=ground\ndatavar 1 groundsize\ntexturevar 0\npolyorivar 2\ntexture -a 1 %s \n0 0 0    1 10  1 0 0    0 1 0\neval poly on\neval polysize %f\neval lum groundsize 0 10\neval polysides 4\neval polymin 0.001 1e10\neval slum 1\neval alpha 1\neval textures on\n\n", $ground, $groundsize/2;  # must divide size by 2


if ($boxheight > 0)
{
    printf CF "object g3=box\n  box -2500,2500 -2500,2500 100,%d\n eval color rgb 0 1 1\n\n", $boxheight;
}

print CF "eval censize 0\neval step 1\neval speed 16\neval run\n\neval censize 0\neval clip 1 1e8\neval jump -25278.4 16496.7 5460.46 -32.5321 -77.908 -102.344 0.999999\n";
print CF "eval focalpoint on\n#eval detach full\neval stereo -0.005 left\n";

#######


sub isnumeric
{    
    my ($x) = @_;    
    my $isnum = 0;    
    if (!($x =~ /[^eE0-9\-\+\.]/))    
    {        
	$isnum = 1;    
    }    
    return $isnum;
}


sub mod 
{
    my($a,$b) = @_;
    return $a % $b;
}

sub min
{
    my($a,$b) = @_;
    return $b if ($a > $b);
    return $a;
}


sub thickline
{
# gets 
#   x0 y0 z0 x1 y1 z1 w c        
# returns string with partiview mesh command
# 
# to join points (x0,y0,z0) and (x1,y1,z1) with a line of thickness w
#  and color c (the index to some cmap file, but this program doesnt
#  have to take care of cmap stuff)

    my ($x0,$y0,$z0,$x1,$y1,$z1,$w,$c) = @_;
    my $answer = "";

    # draw one rectangle

    # first corner is at (a0,b0,c0)

#    $w = 0;

    $a0 = $x0 + $w/2;
    $b0 = $y0 + $w/2;

    if ($z0 == $z1)
    {
	return $answer;
    }

    $c0 = $z0 + $w * ($x1-$x0+$y1-$y0) / (2*($z0-$z1));
    $wcur = sqrt ( ($z0-$c0)**2 + ($y0-$b0)**2 + ($x0-$a0)**2 );

    if ($wcur == 0)
    {
	return $answer;
    }

    #    print "$w $wcur\n";

    $a0 = $x0 + (1.0*$w/$wcur)*($a0-$x0);
    $b0 = $y0 + (1.0*$w/$wcur)*($b0-$y0);
    $c0 = $z0 + (1.0*$w/$wcur)*($c0-$z0);

    # second corner
    
    $a1 = $x0 -($a0-$x0);
    $b1 = $y0 -($b0-$y0);
    $c1 = $z0 -($c0-$z0);

    #third corner

    $a2 = $a1 + ($x1-$x0);
    $b2 = $b1 + ($y1-$y0);
    $c2 = $c1 + ($z1-$z0);

    #fourth corner

    $a3 = $a0 + ($x1-$x0);
    $b3 = $b0 + ($y1-$y0);
    $c3 = $c0 + ($z1-$z0);

    $answer = $answer . "mesh -c $c -s solid {\n2 2\n";
    $answer = $answer . "$a0 $b0 $c0\n";
    $answer = $answer . "$a1 $b1 $c1\n";
    $answer = $answer . "$a3 $b3 $c3\n";
    $answer = $answer . "$a2 $b2 $c2\n";

    $answer = $answer . "}\n";
#    $answer = $answer . "$a0 $b0 $c0 2\n";
#    $answer = $answer . "$a1 $b1 $c1 2\n";
#    $answer = $answer . "$a3 $b3 $c3 2\n";
#    $answer = $answer . "$a2 $b2 $c2 2\n";

    # draw second rectangle

    # first corner is at (p0,q0,r0)
    # (p,q,r)-(x0,y0,z0) is orthogonal to both 
    # (x1,y1,z1)-(x0,y0,z0)  and (a0,b0,c0)-(x0,y0,z0) 

    $p0 = $x0 + ( ($y1-$y0)*($c0-$z0) - ($z1-$z0)*($b0-$y0) );
    $q0 = $y0 + ( ($z1-$z0)*($a0-$x0) - ($x1-$x0)*($c0-$z0) );
    $r0 = $z0 + ( ($x1-$x0)*($b0-$y0) - ($y1-$y0)*($a0-$x0) );
    $wcur = sqrt ( ($z0-$r0)**2 + ($y0-$q0)**2 + ($x0-$p0)**2 );

    if ($wcur == 0)
    {
	return $answer;
    }

#    print "$w $wcur\n";

    $p0 = $x0 + ($w/$wcur)*($p0-$x0);
    $q0 = $y0 + ($w/$wcur)*($q0-$y0);
    $r0 = $z0 + ($w/$wcur)*($r0-$z0);

    # second corner

    $p1 = $x0 -($p0-$x0);
    $q1 = $y0 -($q0-$y0);
    $r1 = $z0 -($r0-$z0);

    #third corner

    $p2 = $p1 + ($x1-$x0);
    $q2 = $q1 + ($y1-$y0);
    $r2 = $r1 + ($z1-$z0);

    #fourth corner

    $p3 = $p0 + ($x1-$x0);
    $q3 = $q0 + ($y1-$y0);
    $r3 = $r0 + ($z1-$z0);

    $answer = $answer . "mesh -c $c -s solid {\n2 2\n";
    $answer = $answer . "$p0 $q0 $r0\n";
    $answer = $answer . "$p1 $q1 $r1\n";
    $answer = $answer . "$p3 $q3 $r3\n";
    $answer = $answer . "$p2 $q2 $r2\n";
    $answer = $answer . "}\n";


#    $answer = $answer . "$p0 $q0 $r0 1\n";
#    $answer = $answer . "$p1 $q1 $r1 1\n";
#    $answer = $answer . "$p3 $q3 $r3 1\n";
#    $answer = $answer . "$p2 $q2 $r2 1\n";

    $L1 = sqrt ( ($a0-$a1)**2 + ($b0-$b1)**2 + ($c0-$c1)**2);
    $L2 = sqrt ( ($p0-$p1)**2 + ($q0-$q1)**2 + ($r0-$r1)**2);
#    print "$L1 $L2\n";

    $L1 = sqrt ( ($a2-$a3)**2 + ($b2-$b3)**2 + ($c2-$c3)**2);
    $L2 = sqrt ( ($p2-$p3)**2 + ($q2-$q3)**2 + ($r2-$r3)**2);
#    print "$L1 $L2\n";
    
#    print $answer,"\n";
    return $answer;
}

sub thinline
{
# gets 
#   x0 y0 z0 x1 y1 z1 w c        
# returns string with partiview mesh command
# 
# to join points (x0,y0,z0) and (x1,y1,z1) with a line of thickness w
#  and color c (the index to some cmap file, but this program doesnt
#  have to take care of cmap stuff)

    my ($x0,$y0,$z0,$x1,$y1,$z1,$w,$c) = @_;
    my $answer = sprintf ("mesh -c %d -s wire {\n1 2\n%0.10f %0.10f %0.10f \n%0.10f %0.10f %0.10f \n}\n", $c,$x0,$y0,$z0,$x1,$y1,$z1);

#    $answer = $answer."\n $x0 $y0 $z0 1 1 $c\n";    # adds a dot at the start of each shower to make it possible to pick (and change the center of rotation) in partiview


    return $answer;
}