#! /usr/bin/perl5

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

$H0 = 75;

while(<>) {
  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(' ', $_);


  next if $isoAr < 10.0;

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

  $ellip{ int(10*$bovera) }++;

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

  $diamkpc = $distmpc * $isoAr * .004848;	# .004848 = 1000 kpc/Mpc / (3600*180/$pi arcsec/radian)
##  printf "%d %.3f %.4g %.4f\n", $type, $bovera, $isoAr, $radius50_iso;
  $name = sprintf("p%d", $type);
  $name =~ tr/-/_/;
  push(@$name, $bovera);
  $var{$type} = $name;
}

print <<EOF;
# isoB_r/isoA_r quantiles as a function of eClass (+/- 0.05 buckets)
# for galaxies of size isoA_r >= 10.0 arcsec.
# Each line has (where q.5 is the median B/A, q.98 the 98th centile etc.):
#   eClass   q.02 q.05 q.25 q.50 q.75 q.95 q.98  total_in_this_eClass

EOF

@pct = (.02, .05, .25, .50, .75, .95, .98);

foreach $type (sort { $a <=> $b } keys(%var)) {
  $name = $var{$type};
  @v = sort( @$name );
  $n = @v;
  printf "eclass %.1f", $type/10;
  for($i = 0; $i < @pct; $i++) {
    printf " %.3g", $v[ int(($n-1) * $pct[$i] + .5) ];
  }
  printf "  %d", $n;
  print "\n";
}

print <<EOF;

# Numbers of galaxies of various isoB/isoA ratios,
# over all galaxies with isoA_r >= 10.0 arcsec.

EOF

foreach $i ( sort { $a<=>$b } keys(%ellip) ) {
  printf "bovera %.1f %d\n", $i/10, $ellip{$i};
}
