#! /usr/bin/perl

&init_eq2gal;

$boxmax = 1e10;
$H0 = 100;
$binstep = .5;

foreach $_ ($0, @ARGV) {
  $cmdline .= sprintf(/\s/ ? ' "%s"' : ' %s', $_);
}

if(@ARGV == 0 && -t STDIN) {
  print STDERR "Usage: $0 [-l labels] [-b rbox(km/s)] [-b xmin,xmax,ymin,ymax,zmin,zmax]
	[-s H0]  nofinger6.short-like-file ... > file.speck
Convert such a file into a .speck list, readable by svd.
With -b, extracts an orthobox:
  -b rbox  selects only galaxies within -rbox <= x,y,z <= +rbox
		of the origin.
  -b xmin,xmax,... selects only those within given limits
		Coordinates given in km/s (i.e., as written in the input file,
		*not* as after scaling to cave units as taken by svd.)
  -s H0	     scale distances and sizes for that value of H0 (default 100).
  -l labels  file of labels, one per line, ``pgcno text text-of-label''
		or ``@ x y z  text  text-of-label''.  Multiple -l options
		read multiple label files.  With no -l, uses built-in defaults.
  -t textures   file of: PGC-or-nNGC-number imagefile
  -d densityfile  sample density function (file written by tully2raw) at each gal
";
  exit(1);
}


while($ARGV[0] =~ /^-[sbhH]/) {
    $nleft = @ARGV;
    shift, $rbox = shift if $ARGV[0] =~ /-b/;
    shift, $H0 = shift if $ARGV[0] =~ /-[shH]/;
    $TEST = shift if $ARGV[0] =~ /^-T/;
    last unless $nleft > @ARGV;
}

@B = grep( $_ ne "", split(/[ ,]+/, $rbox) );
if(@B == 6) {
  ($xmin,$xmax, $ymin,$ymax, $zmin,$zmax) = @B;
} elsif(@B == 1) {
  $xmin = $ymin = $zmin = -$B[0];
  $xmax = $ymax = $zmax = $B[0];
} elsif(@B){
  print STDERR "Expected -b xmin,xmax,ymin,ymax,zmin,zmax OR -b radius; what's \"$rbox\"?\n";
  exit(1);
}


$txno = 30;

foreach $_ (@known) {
  s/\s*#.*//;
  next if $_ =~ /^\s*$/;

  if(/^\s*texturebase\s*(\d+)/) {
    $txno = $1;
    next;
  }

  if(/^\s*\@(.*)/) {
    $_ = $1;
    
    push(@insert, $_); # Anything after an "@", include literally in output.
    next;
  }

  ($pgc, @F) = split(' ', $_);
  if($pgc =~ /^NGC/i) {
    $pgc = "N" . shift(@F) if($pgc =~ /^NGC$/i);
    $pgc =~ tr/a-z/A-Z/;
    $pgc =~ s/NGC/N/;
  } elsif($pgc =~ /^N\d\w*$/i) {
    $pgc =~ tr/a-z/A-Z/;
  }

  next if $pgc eq "";
  while(@F) {
    $key = shift(@F);
    if($key eq "C") {
	$c{$pgc} = shift(@F);
    } elsif($key eq "N") {
	$name{$pgc} = shift(@F);
    } elsif($key eq "P") {
	$pos{$pgc} = shift(@F);
    } elsif($key eq "T" || $key eq "text") {
	($text{$pgc} = join(" ", @F)) =~ tr/_/ /;
	@F = ();
    } elsif($key eq "tx" || $key eq "texture") {
	$texture{$pgc} = $txno;
	push(@insert, sprintf("texture %d %s", $txno++, shift(@F)));
    } elsif($key eq "orient") {
	print STDERR "Expected 6 or 9 numbers (Xx Xy Xz Yx Yy Yz [Px Py Pz]) after \"$pgc $key\"\n"
		unless @F == 6 || @F == 9;
	$orient{$pgc} = join(" ", @F[0..5]);
	$orient{$pgc} .= ";" . join(" ", @F[6..8]) if @F >= 9;
	
    }
  }
}

if($densityfile) {
  open(DENS, $densityfile);
  $_ = <DENS>; # "TULLY # ..."
  $_ = <DENS>; # data...
  @F = split(' ');
  ($dnx, $dny, $dnz) = split(/,/, $F[0]);
  ($drx, $dry, $drz) = split(/,/, $F[1]);
  ($dcx, $dcy, $dcz) = split(/,/, $F[2]);
  ($dscl, $dnommax, $dtruemax) = @F[3..5];
  local($/) = undef;
  $density = <DENS>;
  if($dnz == 0 || $drz == 0 || length($density) < $dnx*$dny*$dnz) {
    print STDERR "$densityfile: Expected tully2raw file\n";
    $density = undef;
  } else {
    $dx0 = $dcx - $drx;  $ddx = 2*$drx/($dnx-1);
    $dy0 = $dcy - $dry;  $ddy = 2*$dry/($dny-1);
    $dz0 = $dcz - $drz;  $ddz = 2*$drz/($dnz-1);
  }
}

$first = 1;

# Map galaxy types to minimum ellipticities -- for estimating inclination
%minba = split(' ', <<EOF);
-6 .5
-5 .5
-4 .45
-3 .4
-2 .2
-1 .2
0  .2
1  .18
2  .18
3  .15
4  .12
5  .1
6  .1
7  .1
8  .1
9  .13
10 .2
EOF

print "#", $cmdline, "\n";


print <<EOF;
# x y z type diam -Mblue txno orient[4..9] cloud+100 # ra+dec name

datavar 0 type
datavar 1 diamkpc
datavar 2 Mg
datavar 3 colorg_r
datavar 4 z
datavar 5 lum
datavar 6 txno
datavar 7 Riso_50
datavar 8 orient0

texturevar txno
polyorivar orient0

EOF

print "\n# Literal labels\n", join("\n", @insert), "\n"
  if(@insert);

while(<>) {
  last if /BEGIN OUTPUT/;
}

while(<>) {
  $patching = ($1 eq "on" || $1 eq "1"), next
	if /^#patch\s*(\S*)/;

  next if /^\s*#|^\s*$/;

# "ra" "dec" "z" "obj.obj.modelMag[0]" "obj.obj.modelMag[1]" "obj.obj.modelMag[2]" "obj.obj.modelMag[3]" "obj.obj.modelMag[4]" "eClass" "obj.petroR50_r" "obj.obj.isoA[2]" "obj.obj.isoB[2]" "obj.obj.isoPhi[2]"

    ($radeg, $decdeg, $redshift, $mu, $mg, $mr, $mi, $mz,
	$eclass, $petroR50r, $isoAr, $isoBr, $isoPhir
    ) = split(' ', $_);


  $distmpc = $redshift * 3.0e5 / $H0;	# Yeah, I know that's not right.
  ($x, $y, $z) = &eq2vec( $radeg, $decdeg );

  $x0 = $x if $x0 > $x;
  $y0 = $y if $y0 > $y;
  $z0 = $z if $z0 > $z;
  $x1 = $x if $x1 < $x;
  $y1 = $y if $y1 < $y;
  $z1 = $z if $z1 < $z;

  $xyz = sprintf("%.5g %.5g %.5g", $x*$distmpc, $y*$distmpc, $z*$distmpc);

  $rahms = &rad2dms( $radeg*$pi/180 / 15 );
  $decdms = &rad2dms( $decdeg*$pi/180 );

  $modulus = log($distmpc/10e-6) * 2.1715;	# * 5/log(10);
  $lumg = exp($modulus - $mg - 19) * .921;	# * log(100)/5

  $type = int($eclass*10 + 10) - 10;	
  $pa = $isoPhir;
  $bovera = $isoBr / $isoAr;

  $radius50_iso = ($petroR50r!=0 ? $isoAr / $petroR50r : 1);

  $diamkpc = $distmpc * $isoAr * .004848;	# .004848 = 1000 kpc/Mpc / (3600*180/$pi arcsec/radian)

    if($type < 0 || $type >= 10) {
	# elliptical galaxies always face the viewer.
	$galeq = "9 9 9 9 9 9";

    } elsif($pa =~ /\d/ && $bovera > 0) {
	@eqbasis = &radec2eqbasis( $radeg, $decdeg );

	# Cook up an inclination from the b/a apparent ellipticity.
	# For the %minba table above, we plotted galaxy-type vs. inclination
	# and, for each type, chose an ellipticity value that most galaxies exceeded.
	# That's assumed to be the edge-on apparent ellipticity.
	# There's probably some clever way of estimating true inclination from
	# these numbers, but here I'll just subtract:
	#   cos i = b/a - { min b/a for this galaxy type }.

	$cosi = $bovera - $minba{ sprintf("%d",$type) };
	$cosi = 0 if $cosi < 0;
	$cosi = 1 if $cosi > 1;
	## $cosi = -$cosi if int($pgc) % 2 == 1;
	$sini = sqrt(1 - $cosi*$cosi);

	$cospa = cos($pa*($pi/180));
	$sinpa = sin($pa*($pi/180));
	# This is really the product of 4 matrices:
	#  (X-rotation by inclination) * (Z-rotation by PA)
	#    * (sky plane orientation)
	# @inclpa is just the product of the first two.
	@inclpa = (
		$cospa,		-$sinpa,	0,
		$cosi*$sinpa,	$cosi*$cospa,	$sini,
		-$sini*$sinpa,	-$sini*$cospa,  $cosi
	    );
	@galeq = &m3mmul( @inclpa, @eqbasis );
	$galeq = sprintf("%.2f %.2f %.2f  %.2f %.2f %.2f", @galeq[0..5]);
	
    } else {
	$galeq = "1 0 0  0 1 0";
    }

    $decdms =~ tr/ /+/;

    $stuff = sprintf("%s %d %.4g %.1f %.2f %.4f %.3g %d %.2f %s # %s%s%s%s",
	$xyz, $type, $diamkpc, $modulus - $mg, $mg - $mr, $redshift, $lumg,
	($type+10), $radius50_iso, $galeq,
	substr($rahms,1,2), substr($rahms,4,2),
	substr($decdms,0,3), substr($decdms,4,2));

    $stuff =~ tr/ / /s;
    print $stuff, "\n";
    $specks++;
}

printf "# Total %d gals, x %.5g %.5g  y %.5g %.5g  z %.5g %.5g (in Mpc, assuming H0 ~ %g)\n\n",
	$specks,
	$x0,$x1,
	$y0,$y1,
	$z0,$z1,
	$H0;



print "\n\n";

print "\n";

while(@records) {
  print shift(@records), "\n";
}


sub dms2rad {
  return ($_[0] + ($_[1] + $_[2]/60)/60) * $pi/180;
}

sub rad2dms {
  local($d) = @_;
  local($sign) = $_[1] || " ";
  $sign = "-", $d = -$d if $d<0;
  $d *= 180/$pi;
  return sprintf("%s%02d:%02d", $sign, int($d), 60*($d - int($d)));
}

sub eq2vec {
  local($ra, $dec) = @_;	# Both RA and DEC in degrees!
  $ra *= $pi/180;
  $dec *= $pi/180;
  local($cdec) = cos($dec);
  return ( cos($ra)*$cdec, sin($ra)*$cdec, sin($dec) );
}

sub cross {
  return ($_[1]*$_[5] - $_[2]*$_[4],
	  $_[2]*$_[3] - $_[0]*$_[5],
	  $_[0]*$_[4] - $_[1]*$_[3]);
}
sub dot {
  $_[0]*$_[3] + $_[1]*$_[4] + $_[2]*$_[5];
}
sub normalize {
   local(@v) = @_;
   local($r) = sqrt(&dot(@v,@v));
   $r=1, $v[0] = 1 if $r == 0;
   return ($v[0]/$r, $v[1]/$r, $v[2]/$r) if @v == 3;
   grep(($_ /= $r) || 1, @v);
}

sub sg2eq {
  &eq2dms( &vm3mul( @_, @Tse ) );
}

sub eq2sg {
  &vm3mul( @_, @Tes );
}

sub radec2eqbasis {
  local($ra, $dec) = @_;
  $ra /= 57.2958;
  $dec /= 57.2958;
  		# vector to galaxy (equatorial coords)
  local(@z) = (cos($ra)*cos($dec), sin($ra)*cos($dec), sin($dec));
  # project @z out of (0,0,1)
  local(@x) = &normalize( -$z[0]*$z[2], -$z[1]*$z[2], 1 - $z[2]*$z[2] );
  local(@y) = &cross( @z, @x );
  (@x, @y, @z);
}

sub eq2dms {
  local($r, $rxy);
  local(@eq) = @_;
  local($ra, $dec);
  $rxy = sqrt($eq[0]*$eq[0] + $eq[1]*$eq[1]);
  $dec = atan2($eq[2], $rxy);
  $ra = atan2($eq[1], $eq[0]);
  $ra += 2*$pi if($ra < 0);
  return(&rad2dms($ra/15), &rad2dms($dec,"+"), sqrt(&dot(@eq,@eq)));
}

sub transpose3 {
  (@_[0,3,6,
      1,4,7,
      2,5,8]);
}

# 3x3 matrix multiply: &mmul(@a, @b) returns @a * @b
sub m3mmul {
  local($i,$j,$k);
  local(@a) = @_[0..8];
  local(@b) = @_[9..17];
  local(@c);
  for($i = 0; $i < 9; $i+=3) {
    $c[$i  ] = $a[$i]*$b[0] + $a[$i+1]*$b[3] + $a[$i+2]*$b[6];
    $c[$i+1] = $a[$i]*$b[1] + $a[$i+1]*$b[4] + $a[$i+2]*$b[7];
    $c[$i+2] = $a[$i]*$b[2] + $a[$i+1]*$b[5] + $a[$i+2]*$b[8];
  }
  return @c;
}

# left-multiply 3-row-vector a by 3x3 matrix T: a*T
sub vm3mul {
  local(@a) = splice(@_,0,3);
  local($i,@v);
  return (
	$a[0]*$_[0] + $a[1]*$_[3] + $a[2]*$_[6],
	$a[0]*$_[1] + $a[1]*$_[4] + $a[2]*$_[7],
	$a[0]*$_[2] + $a[1]*$_[5] + $a[2]*$_[8]);
}


sub init_eq2gal {

  $pi = 3.14159265358979;

  # Both the following matrices are taken from SLALIB routines:
  # 
  # Each row of this matrix is a direction vector
  # in the (L2,B2) galactic system, expressed in J2000 equatorial coordinates.
  #
  @Tge = (
	-0.054875539726,-0.873437108010,-0.483834985808,
	+0.494109453312,-0.444829589425,+0.746982251810,
	-0.867666135858,-0.198076386122,+0.455983795705);
  # 
  # Each column of this matrix is a direction vector in the
  # supergalactic system, expressed in (L2,B2) galactic coordinates.
  #
  @Tsg = (
	-0.735742574804,+0.677261296414,+0.000000000000,
	-0.074553778365,-0.080991471307,+0.993922590400,
	+0.673145302109,+0.731271165817,+0.110081262225);

  @Tse = &m3mmul(@Tsg,@Tge);
  @Tes = &transpose3(@Tse);
  @Teg = &transpose3(@Tge);
}

