Profile: /wise/base/deliv/dev/lib/perl/WISE/FITSIO.pm

File	/wise/base/deliv/dev/lib/perl/WISE/FITSIO.pm
Statements Executed	40
Total Time	0.00972 seconds

Subroutines — ordered by exclusive time
Calls	P	F	Exclusive Time	Inclusive Time	Subroutine
0	0	0	0	0	FITSIO::reform::	area_beams
0	0	0	0	0	FITSIO::reform::	bandnum
0	0	0	0	0	FITSIO::reform::	bandstr
0	0	0	0	0	FITSIO::reform::	ecl2equ
0	0	0	0	0	FITSIO::reform::	ecl2gal
0	0	0	0	0	FITSIO::reform::	inrownum
0	0	0	0	0	FITSIO::reform::	poserr_deg
0	0	0	0	0	FITSIO::reform::	rownum
0	0	0	0	0	FITSIO::reform::	sex2flux
0	0	0	0	0	WISE::FITSIO::	BEGIN
0	0	0	0	0	WISE::FITSIO::	_parse_reform_spec_file
0	0	0	0	0	WISE::FITSIO::	accumhdr
0	0	0	0	0	WISE::FITSIO::	band_merge_simple
0	0	0	0	0	WISE::FITSIO::	find_spatial_cohorts
0	0	0	0	0	WISE::FITSIO::	fits_cat
0	0	0	0	0	WISE::FITSIO::	new
0	0	0	0	0	WISE::FITSIO::	newhdu
0	0	0	0	0	WISE::FITSIO::	sex_to_std
0	0	0	0	0	WISE::FITSIO::	write_bandmerge_file

Line	Stmts.	Exclusive Time	Avg.	Code
1				#! /usr/bin/env perl
2
3				#BEGIN { $Exporter::Verbose=1 }
4
5				package WISE::FITSIO;
6
7	3	5.0e-5	1.7e-5	use WISE::Env (cfglib=>'/wise/base/deliv/dev/lib/perl',import=>[qw/$basedir/]); # spent 593µs making 1 call to WISE::Env::import, max recursion depth 1
8
9	3	2.6e-5	8.7e-6	use strict; # spent 8µs making 1 call to strict::import
10	3	3.1e-5	1.0e-5	use warnings; # spent 24µs making 1 call to warnings::import
11
12	3	3.3e-5	1.1e-5	use vars qw(@ISA $VERSION @EXPORT @EXPORT_OK %EXPORT_TAGS); # spent 72µs making 1 call to vars::import
13
14	3	5.0e-5	1.7e-5	use Carp qw(croak carp confess cluck); # spent 54µs making 1 call to Exporter::import
15
16	1	1.0e-6	1.0e-6	$VERSION = 1.00;
17
18	3	3.1e-5	1.0e-5	use Exporter; # spent 30µs making 1 call to Exporter::import
19	3	0.00057	0.00019	use FITSIO qw(/.*/); # Import everything # spent 295µs making 1 call to Exporter::import
20
21	3	3.0e-5	1.0e-5	use WISE::Release; # spent 5µs making 1 call to import
22	3	0.00073	0.00024	use WISE::Dumper; # spent 49µs making 1 call to Exporter::import
23
24	3	0.00813	0.00271	use File::Basename; # spent 61µs making 1 call to Exporter::import
25
26	1	1.0e-6	1.0e-6	my $version = '$Id: FITSIO.pm 5619 2009-08-21 20:30:50Z tim $ ';
27
28	1	1.1e-5	1.1e-5	@ISA = qw(FITSIO Exporter);
29
30	1	1.0e-6	1.0e-6	@EXPORT = @FITSIO::EXPORT;
31	1	4.0e-6	4.0e-6	@EXPORT_OK = @FITSIO::EXPORT_OK;
32	1	5.0e-6	5.0e-6	%EXPORT_TAGS = %FITSIO::EXPORT_TAGS;
33				#
34				#######
35				#
36				# Override constructor and HDU initialization to add a release header card
37
38				sub new {
39				my $class = shift;
40				my $this = $class->SUPER::new(@_) or return;
41				my $hdunum = $this->currhdu();
42
43				if($this->writeable() &&
44				($this->numhdus()\|\|0) > 0 &&
45				! $this->key("WRELEASE",{hdunum=>1}) ) {
46				$this->key([[
47				WRELEASE=>(WISE::Release->new()->release() // 'NONE').
48				" / WSDS s/w release tag"
49				]],
50				{hdunum=>1});
51				$this->hdunum($hdunum) if $hdunum;;
52				}
53
54				return $this;
55				}
56
57				sub newhdu {
58				my $this = shift;
59				my $rc = $this->SUPER::newhdu(@_) or return;
60				my $hdunum = $this->currhdu();
61
62				if($this->writeable() &&
63				($this->numhdus()\|\|0) > 0 &&
64				! $this->key("WRELEASE",{hdunum=>1}) ) {
65				$this->key([[
66				WRELEASE=>(WISE::Release->new()->release() // 'NONE').
67				" / WSDS s/w release tag"
68				]],
69				{hdunum=>1});
70				$this->hdunum($hdunum) if $hdunum;;
71				}
72
73				return $rc;
74				}
75
76				#
77				#######
78				#
79				# Do the canonical transformatiom from a sextractor output catalog to
80				# the WISE standard catalog format
81				sub sex_to_std {
82				my $this = shift;
83				my $sxcat = shift;
84				my $stdcat = shift;
85				my $opts = shift \|\| {};
86				my $specfile= $opts->{specfile};
87				$specfile = "$basedir/params/cat_reform_specs/sex_to_std.txt"
88				if ! defined $specfile;
89				my %spec;
90				if($specfile) {
91				print "Opening spec file '$specfile' ...\n" if $opts->{verbose};
92				%spec = _parse_reform_spec_file($specfile,$opts)
93				}
94				if(! $stdcat) {
95				$stdcat = WISE::Utils::pathcomp($sxcat,
96				{type=>$opts->{stdtype} \|\|
97				'stdcat'});
98				}
99				print "Reforming '$sxcat' to '$stdcat' ...\n" if $opts->{verbose};
100				my $fileband = ( WISE::Utils::pathdecomp($sxcat) \|\| {} )->{band};
101				my %info = (
102				band => WISE::Utils::bandnum($opts->{band} \|\| $fileband),
103				);
104				my $stdfits = WISE::FITSIO->reform_tbl($sxcat,$stdcat,
105				{%spec,
106				info=>\%info,
107				%$opts});
108				return $stdfits;
109				}
110
111				sub _parse_reform_spec_file {
112				my $specfile = shift;
113				my $opts = shift \|\| {};
114				my $verbose = $opts->{verbose};
115				my $err = "*** $0/PARSESPEC";
116				my %spec;
117
118				open(my $specfh, "<$specfile")
119				or die "$err: Cannot open spec file '$specfile; $!.\n";
120				my $spectxt;
121				$spectxt = join("", grep { ! /^\s*[#;]/ } <$specfh>);
122				# Parse remove-column list
123				if($spectxt =~ /^ \s* remove \s* \{/smx) {
124				my ($rm) = $spectxt =~ /^ \s* remove \s* \{ \s* ( [^}]+ ) \s* \}/smx;
125				die "$err: Bad spec file format in remove definition.\n"
126				if ! defined $rm \|\| $rm =~ /[=<>`'"]/;
127				$spec{rm} = [ split /[\s,]+/, $rm ] if $rm;
128				}
129				# Parse new=old rename-column list
130				if($spectxt =~ /^ \s* rename \s* \{/smx) {
131				my ($mv) = $spectxt =~/^ \s* rename \s* \{
132				\s* ( (?: [^=}\s,]+ \s* = \s* [^=}\s,]+ \s* ,
133				\s* )*
134				(?: [^=}\s,]+ \s* = \s* [^=}\s,]+ )?
135				) \s*
136				\}/smx;
137				die "$err: Bad spec file format in rename definition.\n"
138				if ! defined $mv \|\| $mv =~ /[<>`'"]/;
139				my @mv = split /\s,\s/, $mv if $mv;
140				my %mv;
141				while(@mv) {
142				my($k,$v) = split /\s=\s/, shift(@mv);
143				$mv{$k} = $v;
144				}
145				$spec{mv} = \%mv;
146				}
147				# Parse new=expression or new1,new2,new3=expression add-column list
148				if($spectxt =~ /^ \s* add \s* \{/smx) {
149				my ($add)= $spectxt =~/^ \s* add \s* \{
150				\s* ( \s* [^}]* \s* ) \s*
151				\}/smx;
152
153				die "$err: Bad spec file format in add definition.\n"
154				if ! defined $add;
155				my @add;
156				if($add) {
157				while($add =~ /\s* ([^=}\s]+ \s* = \s* [^=}\s]+) \s* /smxg) {
158				my($k,$v,$o) = split /\s=\s/, $1, 3;
159				die "$err: Bad spec file format at '$k' in add definition.\n"
160				if ! $k \|\| ! defined $v \|\| ! length $v \|\| defined $o;
161				$v =~ s/,\s*$//;
162				push @add, $k, $v;
163				}
164				}
165				$spec{add} = \@add;
166				}
167
168				#print Dumper \%spec;
169
170				return %spec;
171				}
172
173				#
174				#######
175				#
176				# Concatenate FITS tables into a new file
177
178				sub fits_cat {
179				check_sig('$$;$',\@_) or confess;
180				my ($this,$opts,@args) = FITSIO::get_std_args(@_);
181				my $newfits = shift @args;
182				my $fitslist= shift @args;
183				my $like = $opts->{like} \|\| @{ $fitslist \|\| [] }[0];
184				my $cols = $opts->{cols} \|\| '*';
185				my $hdunum = $opts->{hdunum} \|\| 2;
186				my $select = $opts->{select} \|\| '';
187				my $dupras = $opts->{dupr} \|\| 5;
188				my $loncol = $opts->{loncol} \|\| 'ra';
189				my $latcol = $opts->{latcol} \|\| 'dec';
190				my $xcol = $opts->{xcol} \|\| 'x';
191				my $ycol = $opts->{ycol} \|\| 'y';
192				my $rownums = $opts->{rownums};
193				my $rows_in = $opts->{rows};
194				my $addinfo = $opts->{addinfo};
195				my $verbose = $opts->{verbose}\|\| 0;
196				my $nsides = $opts->{nsides};
197				my $nodups = $opts->{nodups};
198				my $extant = $opts->{extant};
199				my $sameisdup=$opts->{sameisdup};
200				my $saverows= $opts->{saverows};
201				my @tracepos= @{ $opts->{trace} \|\| [] };
202				my $dupr = $dupras/3600;
203
204				$\| = 1;
205
206				$rows_in = undef if $rows_in && ! @$rows_in;
207				$rownums = undef if $rownums && ! @$rownums;
208				$nodups = 1 if ! $newfits;
209
210				confess "*** $0/FITSCAT: FITS file list is not an array ref.\n"
211				if $fitslist && ref($fitslist) !~ /array/i;
212
213				print "Concatenating HDU #$hdunum of ".@$fitslist." FITS tables ...\n"
214				if $verbose;
215
216				if($nodups) {
217				my $side = ($dupr2 >= 30 ? $dupr2 : 30) / 3600;
218				print " Cohort resolution initialization for ".
219				"radius = $dupras\" (side=$side deg) ...\n"
220				if $verbose;
221				require WISE::CHealPix;
222				WISE::CHealPix->import('neighbors_ring','side2nsides_deg',
223				'nsides2side_deg', 'ang2pix_ring');
224				require WISE::CoUtils;
225				WISE::CoUtils->import('angdist');
226				$nsides \|\|= side2nsides_deg($side);
227				print " ... HealPix size is ",nsides2side_deg($nsides)*60,
228				"' for nsides=$nsides.\n"
229				if $verbose;
230				}
231
232				confess "*** $0/FITSCAT: No model FITS file defined.\n"
233				if $newfits && ! $extant && ! $like;
234
235				print " Opening FITS output file '$newfits' ...\n"
236				if $verbose && $newfits;
237
238				if(($newfits \|\| $addinfo) && $like) {
239				if( ! ref $like) {
240				my $fits = WISE::FITSIO->new($like) or return;
241				$like = $fits;
242				}
243				} else {
244				# Don't need model
245				$like = undef;
246				}
247
248				if($like && $addinfo) {
249				# Add columns to hold input file info
250				my $model = $like->file();
251				print " Re-opening model '$model' as virtual ".
252				"to add info columns ...\n"
253				if $verbose;
254				my $addcols = "";
255				$addcols = "[".($hdunum-1)."]".
256				"[col FILENUM(J)=0 ; FILENAME(80A)='none' ; *]";
257				$like->end;
258				$like = WISE::FITSIO->new($model.$addcols,{lc=>1})
259				or return;
260				}
261
262				my ($cen,$width,$no_rep_select);
263				if($nodups && $like) {
264				$cen = $like->key("CRPIX1",{hdunum=>1}) or return;
265				$width = ($cen - 0.5)*2;
266				print " Found pixel center,width ($cen,$width).\n"
267				if $verbose > 1;
268				} else {
269				# Since we can't get image details, we must not care about
270				# selecting the best representative from the dups
271				$no_rep_select = 1;
272				}
273
274				my ($fitsout,$atrow);
275				if($newfits) {
276				$fitsout = WISE::FITSIO->new($newfits,
277				{
278				mode=>$extant?'rw':'new',
279				like=>$extant?undef:$like,
280				})
281				or return;
282
283				$fitsout->hdunum($hdunum);
284				$atrow = $extant ? $fitsout->size() + 1 : 1;
285				print " Appending at record $atrow ...\n"
286				if $verbose && $atrow>1;
287				} else {
288				$atrow = 1;
289				}
290
291				return $fitsout if ! $fitslist \|\| ! @$fitslist; # No input; stop here
292
293
294				if($select) {
295				$select = "[".($hdunum-1)."][$select]";
296				print " Using select text '$select'.\n" if $verbose;
297				}
298
299				my $nrowsread = 0;
300				my $nrowsproc = 0;
301				my $filenum = 1;
302				my (@saved_rows,@saved_fits);
303				my (@saved_refs,%pixs); # For dup removal
304				print " Beginning file I/O ...\n" if $verbose;
305				FILE: for my $file (@$fitslist) {
306				my $fits;
307				my $infile;
308				if(! ref $file) {
309				print " Opening input FITS file '$file' (#$filenum) ...\n"
310				if $verbose > 2;
311				$fits = WISE::FITSIO->new($file.$select,{lc=>1}) or return;
312				} else {
313				# Already opened
314				$fits = $file;
315				}
316				$infile = $fits->file();
317				$fits->hdunum($hdunum);
318				my $nrows = $fits->size();
319				my $rows;
320				if(! $rows_in) {
321				print " Reading input FITS file '$file' (#$filenum) ...\n"
322				if $verbose > 1;
323				$rows = $fits->readcols($cols);
324				print " ... read $nrows rows ...\n"
325				if $verbose > 1;
326				return if ! defined $rows;
327				next FILE if $nrows == 0;
328				if($addinfo) {
329				#print join(",",keys%$rows),"\n";
330				$rows->{FILENUM} = [];
331				$rows->{FILENUM}[$_] = $filenum for 0..$nrows-1;
332				my $filebase = basename($infile);
333				$rows->{FILENAME}= [];
334				$rows->{FILENAME}[$_] = $filebase for 0..$nrows-1;
335				}
336				} else {
337				$rows = $rows_in->[$filenum-1];
338				#print "--- ".@{$rows->{$loncol}}.", ".
339				# join(",",keys%{$rownums->[$filenum-1]})."\n";
340				}
341				confess "*** $0/FITSCAT: Table missing one or more positional ".
342				"rows $loncol $latcol, $xcol, $ycol.\n"
343				if grep { ! $rows->{$_} } ($loncol, $latcol, $xcol, $ycol);
344				push @saved_rows, $rows if $saverows;
345				push @saved_fits, $fits if $saverows;
346				if($rownums && ! $rownums->[$filenum-1]) {
347				print" Skipping file $filenum since no rows were requested.\n"
348				if $verbose;
349				next FILE ;
350				}
351				if($fitsout) {
352				$fitsout->writecols($rows,{row=>!$nrowsread?$atrow:undef})
353				or return;
354				}
355				if($nodups) {
356				ROWS: for my $i (0..$nrows-1) {
357				next ROWS if $rownums && $rownums->[$filenum-1] &&
358				! $rownums->[$filenum-1]{$i+1};
359				++$nrowsproc;
360				my $outrow = $atrow + $nrowsread + $i;
361				my $pix = ang2pix_ring($nsides,
362				$rows->{$loncol}[$i],
363				$rows->{$latcol}[$i]);
364				my $inland;
365				if(! $no_rep_select) {
366				my $xin = ($rows->{$xcol}[$i] < $cen
367				? $rows->{$xcol}[$i] - 0.5
368				: $width - $rows->{$xcol}[$i] + 0.5);
369				my $yin = ($rows->{$ycol}[$i] < $cen
370				? $rows->{$ycol}[$i] - 0.5
371				: $width - $rows->{$ycol}[$i] + 0.5);
372				$inland = $xin < $yin ? $xin : $yin;
373				}
374				my $ref = {lon => $rows->{$loncol}[$i],
375				lat => $rows->{$latcol}[$i],
376				inland => $inland,
377				rownum => $outrow,
378				filerow => $i+1,
379				filenum => $filenum,
380				pix => $pix,
381				dups => undef};
382				# Save info for stepping through sequentially later
383				push @saved_refs, $ref;
384				# Index for later retrieval in dup look-up
385				my @pixs = ($pix,neighbors_ring($nsides,$pix));
386				print "--- Pixels for $outrow : ".
387				"$rows->{$loncol}[$i],$rows->{$latcol}[$i] : @pixs\n"
388				if @tracepos &&
389				angdist($rows->{$loncol}[$i],$rows->{$latcol}[$i],
390				@tracepos)*3600<$dupras;
391				push @{ $pixs{$_}\|\|=[] }, $ref for @pixs;
392				} # ROWS
393				} # $noups
394				$nrowsread += $nrows;
395				print " ... through with file $filenum; ".
396				"$nrowsread rows read and processed $nrowsproc so far.\n"
397				if $verbose>1;
398				} continue {
399				++$filenum;
400				} # FILE
401
402
403				print "Concatenation complete. Read $nrowsread rows and ".
404				"processed $nrowsproc from ".@$fitslist." files.\n".
405				($atrow>1
406				?"Output file now has ".($nrowsread+$atrow-1)." rows.\n":"")
407				if $verbose;
408
409				if($nodups) {
410				print " Looking for dups ...\n" if $verbose;
411				for my $ref (@saved_refs) {
412				# Look up spatial index
413				my @cands = @{ $pixs{$ref->{pix}} \|\| [] };
414				my @dups;
415				my $trace = @tracepos &&
416				grep {angdist(@{$_}{'lon','lat'},@tracepos)*3600
417				<$dupras}
418				@cands;
419				print "--- Candidates of $ref->{rownum} @{$ref}{'lon','lat'}: ".
420				join(",",map{"$_->{rownum}"}@cands)."\n"
421				if $trace;
422				# Step through candidates
423				for my $cand (@cands) {
424				# Ignore self
425				next if $ref->{lon}==$cand->{lon} && $ref->{lat}==$cand->{lat};
426				# Match to radius
427				my $r = angdist(@{$ref}{'lon','lat'},
428				@{$cand}{'lon','lat'});
429				print "--- $cand->{rownum}: @{$cand}{'lon','lat'} = ".
430				"$r <= $dupr\n"
431				if $trace;
432				next if $r > $dupr;
433				# Found a dup
434				print "--- ADD $cand->{rownum}\n" if $trace;
435				push @dups, $cand;
436				}
437				$ref->{dups} = \@dups;
438				print "--- Dups: ". join(",",map{"$_->{rownum}"}@dups)."\n"
439				if $trace;
440				}
441				print " Creating dup cohorts ...\n" if $verbose;
442				for my $ref (@saved_refs) {
443				next if $ref->{alldups}; # Already done
444				next if ! @{$ref->{dups}}; # Nothing to do
445				print "--- $ref->{rownum} @{$ref}{'lon','lat'}: ".
446				join(",",map{"$_->{rownum}/$_->{inland}"}
447				@{$ref->{dups}})."\n"
448				if @tracepos &&
449				angdist(@{$ref}{'lon','lat'},@tracepos)*3600<$dupras;
450				my (%toadd,%alldups);
451				@toadd{@{$ref->{dups}}} = @{$ref->{dups}};
452				# Iteratively add the dup lists from all dups, and dups of dups ...
453				COHORTS: while(1) {
454				my @toadd = grep {$_} values %toadd;
455				#print "--- @toadd\n";
456				last COHORTS if ! @toadd;
457				for my $cand (@toadd) {
458				$alldups{$cand} = $cand;
459				my @add = grep {! exists $toadd{$_}} @{ $cand->{dups} };
460				@toadd{@add} = @add;
461				$toadd{$cand} = 0; # Mark as done
462				}
463				}
464				my @alldups = values %alldups;
465				$_->{alldups} = \@alldups for @alldups;
466				print "--- $ref->{rownum} @{$ref}{'lon','lat'}: ".
467				join(",",map{"$_->{rownum}/$_->{inland}"}
468				@alldups)."\n"
469				if @tracepos &&
470				angdist(@{$ref}{'lon','lat'},@tracepos)*3600<$dupras;
471				}
472				my @duprows;
473				if(! $no_rep_select) {
474				print " Selecting best dup from cohorts ...\n" if $verbose;
475				my %duprows;
476				for my $ref (@saved_refs) {
477				next if ! $ref->{alldups}; # No dups, or cohort already done
478				my $maxref;
479				my @alldups = @{$ref->{alldups}};
480				print "--- Cohort:\n\t".
481				join("\n\t",map{"$_->{rownum}:$_->{lon},$_->{lat}"}
482				@alldups).
483				"\n"
484				if @tracepos &&
485				angdist(@{$ref}{'lon','lat'},@tracepos)*3600<$dupras;
486				for my $cand (@alldups) {
487				$maxref = ! $maxref \|\| $cand->{inland}>$maxref->{inland}
488				? $cand : $maxref;
489				delete $cand->{alldups}; # Remove so it's processed once
490				$duprows{$cand->{rownum}} = $cand->{rownum}; # Mark for rm
491				print "--- removing $cand->{rownum} ($cand->{inland})\n"
492				if @tracepos &&
493				angdist(@{$cand}{'lon','lat'},@tracepos)
494				*3600<$dupras;
495				}
496				# Allow winner to survive
497				delete $duprows{$maxref->{rownum}};
498				print "--- resurrecting $maxref->{rownum}".
499				#"; ".join(",",sort {$a<=>$b} keys %duprows).
500				"\n"
501				if @tracepos &&
502				angdist(@{$ref}{'lon','lat'},@tracepos)*3600<$dupras;
503				if(! $sameisdup) {
504				# Also rescue those on the same image
505				delete $duprows{$_->{rownum}}
506				for grep {exists $duprows{$_->{rownum}} &&
507				$maxref->{filenum}==$_->{filenum}} @alldups;
508				}
509				}
510				@duprows = sort {$a<=>$b} keys %duprows;
511				print " Found ".@duprows." duplicate rows.\n"
512				if $verbose;
513				} # rep selection
514				if(@duprows && $fitsout) {
515				print " Deleting duplicate rows ...\n"
516				if $verbose;
517				# Close and reopen
518				$fitsout->end();
519				$fitsout = WISE::FITSIO->new($newfits, {mode=>'rw'})
520				or return;
521				@duprows = sort {$a <=> $b} @duprows;
522				$fitsout->rmrows(\@duprows,{hdunum=>$hdunum})
523				or return;
524				}
525				} # nodups
526
527				if($fitsout) {
528				print "Final output file has ".$fitsout->size()." rows.\n"
529				if $verbose;
530				return wantarray ? ($fitsout, \@saved_refs, \@saved_rows, @saved_fits,
531				\%pixs, $nsides)
532				: $fitsout;
533				} else {
534				# No output file requested, so just return the data,
535				# including dup cohort data and index
536				return wantarray ? (\@saved_refs, \@saved_rows, \@saved_fits,
537				\%pixs, $nsides)
538				: \@saved_refs;
539				}
540
541				}
542
543				# Front-end for above that just wants the dup cohort list
544
545				sub find_spatial_cohorts {
546				check_sig('$;$',\@_) or confess;
547				my ($this,$opts,@args) = FITSIO::get_std_args(@_);
548				my $fitslist= shift @args;
549				$opts->{dupr} = $opts->{r} if $opts && $opts->{r};
550				$fitslist = [ $fitslist ] if ! ref $fitslist;
551				return $this->fits_cat(0,$fitslist,$opts);
552				}
553
554				sub band_merge_simple {
555				check_sig('$;$',\@_) or confess;
556				my ($this,$opts,@args) = FITSIO::get_std_args(@_);
557				my $fitslist= shift @args;
558				$opts->{loncol} \|\|= 'ra';
559				$opts->{latcol} \|\|= 'dec';
560				$opts->{xcol} \|\|= 'x';
561				$opts->{ycol} \|\|= 'y';
562				print "Simple band merge ...\n" if $opts->{verbose};
563				# Assume files represent bands (one per) we wish to merge
564				# Get all cohorts
565				my $nit = 1;
566				my (@bandnums,@merges, $cohorts, $rows, $fits, @redo_rows, $svrows);
567				push @bandnums, WISE::Utils::bandnum(
568				WISE::Utils::pathdecomp($_)->{band})
569				for @$fitslist;
570				my %bandfile = map { ($bandnums[$_] => $_) } 0..$#bandnums;
571				my %fileband = reverse %bandfile;
572				COHORTS: { # Outer loop to iterate as necessary to process leftovers
573				print " Building cohorts ($nit) ...\n" if $opts->{verbose};
574				($cohorts, $rows, $fits) =
575				$this->find_spatial_cohorts($fitslist,
576				{saverows=>1,
577				rows=>$svrows,
578				rownums=>\@redo_rows,
579				%$opts,
580				verbose=>$opts->{verbose} ? 2 : 0,
581				})
582				or return;
583				@redo_rows = ();
584				my $nreproc = 0;
585				# Step through cohorts and build merge list, breaking ambiguities where
586				# necessary.
587				print " Merging ...\n" if $opts->{verbose};
588				for my $cohort (@$cohorts) {
589				next if $cohort->{done}; # Already in another merge item
590				# Remove superfluous structure element
591				$cohort->{dups}=undef;
592				push(@merges, [$cohort]), next # Trivial merge; no dups
593				if ! $cohort->{alldups};
594				# Trivial merge; self only
595				$cohort->{alldups}=undef, push(@merges, [$cohort]), next
596				if @{$cohort->{alldups}}==0 \|\|
597				(@{$cohort->{alldups}}==1 && $cohort->{alldups}[0]==$cohort);
598				# There is a non-trivial cohort list.
599				$cohort = $cohort->{alldups}; # Save a redirection
600				$_->{dups}=$_->{alldups}=undef for @$cohort; # Avoid self-pointing
601				$_->{done} = 1 for @$cohort; # Mark as done
602				# Classify by file of origin.
603				my %byfile;
604				push @{ $byfile{$_->{filenum}} }, $_ for @$cohort;
605				# If all members come from different files, we're good
606				push(@merges, $cohort), next
607				if grep(@{$byfile{$_->{filenum}}}==1, @$cohort) == keys %byfile;
608				# If only one band is involved, just split them into individual
609				# cohorts and move on
610				if(keys(%byfile) == 1) {
611				$_->{ambiguous} = -1 for @$cohort; # Mark as ambiguous
612				push @merges, [$_] for @$cohort;
613				next;
614				}
615				# We have a complex situation
616				$_->{ambiguous} = -2 for @$cohort; # Mark as even more ambiguous
617				# Break ambiguity for bands/files with multiple matches
618				# Count number of files/bands with multiple cohort members since this
619				# may (eventually) determine how we break up the group
620				my @multfiles = grep { @{ $byfile{$_} } > 1 } keys %byfile;
621				my @rm;
622				# Step through first source from each file with multiple members
623				for my $src (map { $byfile{$_}[0] } @multfiles) {
624				my @infile = @{ $byfile{$src->{filenum}} }; # All srcs in this file
625				# Stupid tie-breaking rule: Get the average distance from
626				# each source in the band to the other sources in the
627				# cohort one with the smallest average distance
628				my ($minaver,$mindup);
629				for my $dup ( @infile ) { # Step through others in same file/band
630				$dup->{ambiguous} = 1+@infile;
631				my $aver = 0;
632				my $n = 0;
633				for my $src2 (@$cohort) {
634				next if $src2 == $dup; # Not self
635				# If it's a real cluster fsck (ALL bands have
636				# multiple members in the cohort), all we can do
637				# (with this simple algorithm) is compare to all
638				# cohort members. But if any files/bands are
639				# singlette, just compare to them.
640				next if @multfiles < keys %byfile && # singlettes exist
641				@{ $byfile{$src2->{filenum}} } > 1; # we're not one
642				$aver += angdist(@{$dup}{'lon','lat'},
643				@{$src2}{'lon','lat'} )*3600;
644				++$n;
645				}
646				$aver /= $n;
647				$dup->{aver} = $aver;
648				if(! defined $minaver \|\| $aver <= $minaver) {
649				$minaver = $aver;
650				$mindup = $dup;
651				}
652				}
653				# Now $mindup is the "best" representative for this band
654				# Break off the other dups into single-band "floaters" and
655				# either a) remember their rows to reprocess them if
656				# multiple bands had multiple cohort members, or b) just
657				# make them into cohorts of their own if not.
658				for my $dup ( @infile ) {
659				next if $mindup == $dup; # Skip the chosen one
660				if(@multfiles < 2) {
661				# Add as a new, single-member cohort
662				$dup->{ambiguous} = -3;
663				push @merges,[$dup];
664				} else {
665				# Breal off and reprocess on a subsequent iteration
666				$redo_rows[$dup->{filenum}-1] \|\|= {};
667				$redo_rows[$dup->{filenum}-1]{$dup->{filerow}+0} = 1;
668				++$nreproc;
669				}
670				push @rm, $dup; # Mark for removal from cohort
671				}
672				}
673				# Remove dups we cut loose from the herd
674				for my $rm (@rm) {
675				# Find index in cohort list
676				my ($ix) = grep {$cohort->[$_] == $rm} 0..$#$cohort;
677				# Remove it from the cohort
678				splice @$cohort,$ix,1;
679				}
680				# Push final resolved cohort
681				push @merges, $cohort;
682				}
683				if($opts->{verbose}) { # Print stats
684				my $nbandsrcs;
685				my %nbands;
686				my %nbandcombos;
687				my %nambig;
688				for my $mrg (@merges) {
689				$nbandsrcs += @$mrg;
690				++$nbands{@$mrg};
691				my $combo = "";
692				$combo = join "", sort map {$fileband{$_->{filenum}-1}} @$mrg;
693				++$nbandcombos{$combo};
694				++$nambig{$_->{ambiguous}\|\|0} for @$mrg;
695				}
696				print " Total mergers = ".@merges." representing $nbandsrcs ".
697				"detections; Will reprocess $nreproc.\n";
698				print " Band count hist: ".join(", ",
699				map {"$_:$nbands{$_}" }
700				sort keys %nbands)."\n";
701				print " Band combo hist: ".join(", ",
702				map {"$_:$nbandcombos{$_}" }
703				sort keys %nbandcombos)."\n";
704				print " Ambiguity hist : ".join(", ",
705				map {"$_:$nambig{$_}" }
706				sort {$a<=>$b} keys %nambig)."\n";
707				}
708				if(@redo_rows) {
709				# Rerun cohort finding on reduced set of rows to resolve sticky wickets
710				print "Reprocessing ".
711				join(", ",map {scalar(keys %{$redo_rows[$_]})." rows from file ".
712				($_+1)}
713				grep {$redo_rows[$_]} 0..$#redo_rows)."\n"
714				if $opts->{verbose};
715				$fitslist = $fits; # Avoid reopening
716				$svrows = $rows; # Avoid rereading
717				++$nit;
718				redo COHORTS;
719				}
720
721				} # COHORT
722				print "Merging complete.\n" if $opts->{verbose};
723				return wantarray ? (\@merges, $rows, $fits) : \@merges;
724				}
725
726				sub write_bandmerge_file {
727				check_sig('$$;$',\@_) or confess;
728				my ($this,$opts,@args) = FITSIO::get_std_args(@_);
729				my $mergefile = shift @args;
730				my $filelist = shift @args;
731				my $mergecols = $opts->{mergecols}
732				\|\| [qw/ra dec flux sigflux mag x y sigx sigy snr src/];
733				my $xycols = $opts->{xycols}
734				\|\| [qw/x y/];
735				my $bands = $opts->{bands};
736				my $bandwcs = $opts->{wcs};
737				my $merges = $opts->{merges};
738				my $rows = $opts->{rowdata};
739				my $fits = $opts->{fitsobjects};
740				my $debug = $opts->{debug};
741				my %xycols = map { (lc $_ => 1) } @$xycols;
742				if(! $bands) {
743				$bands = [];
744				for my $file (@$filelist) {
745				push @$bands, WISE::Utils::bandstr(
746				WISE::Utils::pathdecomp($file)->{band});
747				}
748				}
749				my @bandnums = map { WISE::Utils::bandnum($_) or die; } @$bands;
750				my %bandfile = map { ($bandnums[$_] => $_) } 0..$#bandnums;
751
752				if($bandwcs) {
753				require Astro::WCS::LibWCS;
754				}
755
756				#print "--- /@$bands/ /@bandnums/ /@bandnums/\n";
757				print "Writing columns '@$mergecols' from ".(@$filelist).
758				" band files for bands @bandnums to merge file '$mergefile' ...\n"
759				if $opts->{verbose};
760				if(! $merges) {
761				print " Performing merge ...\n" if $opts->{verbose};
762				($merges,$rows,$fits) = $this->band_merge_simple($filelist,$opts)
763				or return;
764				}
765				warn("*** $0/WRTMERGE: Missing row data or FITS objects.\n"), return
766				if ! $rows \|\| ! $fits;
767				warn("*** $0/WRTMERGE: Incorrect row data or FITS object count compared ".
768				"to supplied file list.\n"), return
769				if @$rows != @$filelist \|\| @$fits != @$filelist;
770				print " Contructing new records for merge file ...\n"
771				if $opts->{verbose};
772				my $colinfo = $fits->[0]->colinfo('*',{hdunum=>2});
773				$colinfo->{lc $_} = $colinfo->{$_} for keys %$colinfo;
774				#print Dumper $colinfo;
775				my $data;
776				my $mrgnum = 0;
777				for my $mrg (@$merges) {
778				my %gotband;
779				# Add some overall band info
780				++$mrgnum;
781				my $nbands = @$mrg;
782				my $combonums = join "", sort map { $bandnums[$_->{filenum}-1] } @$mrg;
783				my $combobits = 0;
784				$combobits \|= $_-1 for sort map { $bandnums[$_->{filenum}-1] } @$mrg;
785				push @{ $data->{"merge_id"} }, $mrgnum;
786				push @{ $data->{"nbands"} }, $nbands;
787				push @{ $data->{"bandbits"} }, $combobits;
788				push @{ $data->{"bandnums"} }, $combonums;
789				my $ambig = 0;
790				abs($_->{ambiguous}\|\|0)>abs($ambig) and $ambig=$_->{ambiguous}
791				for @$mrg;
792				push @{ $data->{"ambig"} }, $ambig;
793				# Get representative position; use lowest detection present
794				my $repmrg = (sort {$a->{filenum} <=> $b->{filenum}} @$mrg)[0];
795				my $repfilenum = $repmrg->{filenum};
796				my $reprows = $rows->[$repfilenum-1];
797				my $rownum = $repmrg->{filerow};
798				my ($ra,$dec) = ($reprows->{ra}[$rownum-1],
799				$reprows->{dec}[$rownum-1]);
800				my (@x, @y);
801				if($bandwcs) {
802				my $off = 0;
803				for my $bix (0..$#$bands) {
804				my $b = $bandnums[$bix];
805				my $wcs = $bandwcs->[$bix];
806				$wcs->wcs2pix($ra,$dec,$x[$b-1],$y[$b-1],$off);
807				}
808				#print "x//y = @x // @y\n";
809				}
810				push @{$data->{"ra"} }, $ra;
811				push @{$data->{"dec"}}, $dec;
812				for my $src (@$mrg) {
813				my $fileix = $src->{filenum}-1;
814				my $filerows= $rows->[$fileix];
815				my $bandnum = $bandnums[$fileix];
816				my $rownum = $src->{filerow};
817				#print "/$fileix/$bandnum/$rownum/\n";
818				$gotband{$bandnum} = $fileix;
819				for my $col (@$mergecols) {
820				$data->{"${col}_$bandnum"} \|\|= [];
821				push @{ $data->{"${col}_$bandnum"} },
822				$filerows->{$col}[$rownum-1];
823				}
824				push @{ $data->{"rownum_$bandnum"} }, $rownum;
825				push @{ $data->{"filenum_$bandnum"} }, $fileix+1;
826				}
827				# Fill other bands
828				for my $bix (0..$#$bands) {
829				my $bandnum = $bandnums[$bix];
830				#print "------- $bandnum /$gotband{$bandnum}/\n";
831				next if defined $gotband{$bandnum};
832				#print "------- adding blanks for $bandnum /@$mergecols/\n";
833				for my $col (@$mergecols) {
834				$data->{"${col}_$bandnum"} \|\|= [];
835				my $blank = $colinfo->{lc $col}{type} eq 'A' ? '' : -9999;
836				my $val;
837				if(! $bandwcs \|\| ! $xycols{lc $col}) {
838				$val = $blank;
839				} else {
840				# X/Y column and we can generate a value because we
841				# have wcs info
842				$val = $col =~ /x/i ? $x[$bandnum-1] : $y[$bandnum-1];
843				}
844				push @{ $data->{"${col}_$bandnum"} }, $val;
845				#print "------- /$bandnum/$col/$colinfo->{lc $col}{type}".
846				# "/$val/".@{$data->{"${col}_$bandnum"}}."/\n";
847				}
848				push @{ $data->{"rownum_$bandnum"} }, -9999;
849				}
850				#print join(",", map { "$_=>/".$data->{$_}[-1]."/" }
851				# (qw(MERGE_ID NBANDS BANDS),
852				# (map {my$b=$_; (map{"${_}_$b"} qw(ROW ELON ELAT FLUX))}
853				# @freqstrs))
854				# ),"\n";
855				}
856				print " Open new merge file '$mergefile' ...\n" if $opts->{verbose};
857				# Set up new HDU def. with columns in band number order
858				my @hdudefs;
859				my $colnum = 0;
860				# New global columns
861				my %newtypes = qw(
862				ambig 1J
863				merge_id 1J nbands 1J bandbits 1J bandnums 10A rownum 1J
864				ra 1D dec 1D filenum 1J
865				);
866				for my $new (qw(merge_id nbands bandbits bandnums ambig ra dec)) {
867				++$colnum;
868				die "*** $0/WRTMERGE: Never heard of column '$new'.\n"
869				if ! $newtypes{lc $new};
870				push @hdudefs, {name => "TTYPE$colnum",
871				value=> $new,
872				type => "C"},
873				{name => "TFORM$colnum",
874				value=> $newtypes{lc $new},
875				type => "C"};
876				}
877				for my $bandnum (sort {$a<=>$b} @bandnums) {
878				my $fileix = $bandfile{$bandnum};
879				my $bandnum = $bandnums[$fileix];
880				# New per-band cols
881				for my $new (qw(rownum filenum)) {
882				++$colnum;
883				die "*** $0/WRTMERGE: Never heard of column '$new' ($bandnum).\n"
884				if ! $newtypes{lc $new};
885				push @hdudefs, {name => "TTYPE$colnum",
886				value=> "${new}_$bandnum",
887				type => "C"},
888				{name => "TFORM$colnum",
889				value=> $newtypes{lc $new},
890				type => "C"};
891				}
892				# Old per band cols
893				for my $col (@$mergecols) {
894				++$colnum;
895				die "*** $0/WRTMERGE: Never heard of column '$col' ($bandnum).\n"
896				if ! $colinfo->{lc$col};
897				push @hdudefs, {name => "TTYPE$colnum",
898				value=> "${col}_$bandnum",
899				type => "C"},
900				{name => "TFORM$colnum",
901				value=> $colinfo->{lc$col}{repeat}.
902				$colinfo->{lc$col}{type},
903				type => "C"};
904				}
905				}
906				print " ... adding $colnum column defs ...\n"
907				if $opts->{verbose};
908				#print Dumper \@hdudefs;
909				my $mrgfits = WISE::FITSIO->new($mergefile,
910				{
911				mode => 'write',
912				hdudefs => [
913				# HDU#1
914				[
915				@{ $fits->[0]->gethdu(1) },
916				map {
917				{name=>'INFILE'.($_+1),
918				value=>$filelist->[$_],
919				type=>'C',
920				}
921				} 0..$#$filelist,
922				],
923				# HDU#2
924				[
925				{name => 'XTENSION',
926				value=> 'BINTABLE',
927				type => 'C'},
928				@hdudefs,
929				],
930				],
931				})
932				or return;
933				print " ... writing new columns.\n"
934				if $opts->{verbose};
935				if($debug) {
936				print "----- ";
937				print $_," => ",WISE::Utils::undefize($data->{$_}[0]),", "
938				for sort keys %$data;
939				print "\n";
940				}
941				$mrgfits->writecols($data,{hdunum=>2}) or return;
942				print "Band merge file written.\n" if $opts->{verbose};
943				return wantarray ? ($merges,$data,$rows,$fits) : $merges;
944				}
945
946				#
947				#######
948				#
949				# Update the first HDU data structure in $hdu0 to accumulate info from
950				# the header info in the FITS file object $fits using the accumulation
951				# rules in $rules.
952				# Form of the rules is a hash of pairs like:
953				# NAME=><instructions>
954				# where NAME is the FITS header card name
955				# and <instructions can be
956				# 1 = use the value from the first header
957				# 2 = use the value from the last header
958				# 3 = add the values
959				# 4 = average the values
960				# XXX@@@ = keep all by adding cards named XXX followed by a number with
961				# as many digits as there are @'s.
962				# [<number>,XXX@@@] writing = do both a number rule and a '@' rule.
963				# $rule can also be a list of strings in the above form OR a macro
964				# like '@1' which says look up macro one in %stdrules.
965
966				sub accumhdr {
967				check_sig('$$$;$',\@_) or confess;
968				my ($this,$opts,@args) = FITSIO::get_std_args(@_);
969				my $n = shift @args; # Sequential file number
970				my $rules = shift(@args) \|\| {}; # Header card update rules
971				my $accum = shift(@args) \|\| {}; # Accumulating header info (also returned)
972				my $debug = $opts->{debug} \|\| 0;
973				my $verbose = $opts->{verbose};
974				my %stdrules = %{ $opts->{stdrules} \|\| {} };
975				my @minkeep = @{ $opts->{keep} \|\| [] };
976				my $hdunum = $opts->{hdunum} \|\| 1;
977				my $blank = $opts->{blank};
978				my $add = $opts->{add};
979				my ($hdu,@blank);
980
981				# Initialize the blank value to something impossible to match if it
982				# wasn't passed by the caller. This makes some logic below easier.
983				$blank = "\0" if ! defined $blank; # NULL should never be in a FITS header.
984				@blank = ref $blank ? (grep {defined $_} @$blank) : ($blank);
985
986				print "Blank(s) = /",join("/",@blank),"/\n" if $debug =~ /h/;
987
988				if($debug =~ /h/) { require Data::Dumper; Data::Dumper->import(); }
989
990				# If $rules is a string, It comes (more or less) straight
991				# from the command line; parse it into a hash.
992				if($rules && ! ref $rules) {
993				# Process which header cards to accumulate.
994				# Most often this will be a default set defined by the %stdrules hash.
995				# Look for '@whatever' macros and expand them.
996				while($rules =~ /(?:^\|,)\s*(@[^,]+)/g) {
997				print "Looking at rules macro $1 ...\n" if $verbose;
998				confess "*** $0/accumhdr: Standard card set '$1' not recognized.\n"
999				if ! $stdrules{$1};
1000				my $sub = ref($stdrules{$1}) ? join(",",@{$stdrules{$1}})
1001				: $stdrules{$1};
1002				substr($rules, $-[1], $+[1] - $-[1]) = $sub;
1003				print "... rules now '$rules'\n" if $debug;
1004				}
1005				# Evaluate the full hash in a safe compartment
1006				my %rules = eval "($rules)";
1007				confess "*** $0/accumhdr: Could not interpret rules parameter:\n".
1008				" rules='$rules'.\n $@"
1009				if $@;
1010				#print "rules =\n",Dumper(\%rules) if $debug && %rules;
1011
1012				$rules = 0 if ! %rules;
1013				$rules = \%rules if %rules;
1014
1015				}
1016
1017				confess "*** $0/accumhdr: No accumulation rules defined.\n"
1018				if ! $rules;
1019
1020				# Read current HDU from file
1021				$hdu = $this->gethdu($hdunum);
1022
1023				# Add in new info
1024				if($add && @$add) {
1025				my @add = @$add;
1026				while(@add) {
1027				push @$hdu,{name=>shift(@add),value=>shift(@add)};
1028				}
1029				}
1030
1031				#print "n=$n, new=".($accum\|\|"").", HDU=$hdu\n"
1032				# if $debug =~ /h/;
1033
1034				#print "New HDU =\n",Dumper($hdu) if $debug =~ /hh/;
1035
1036				# Step through all current FITS HDU #1 header cards.
1037				for my $r (@$hdu) {
1038				my $k = $r->{name}; # Get the key name
1039				my $newname = $k; # Changed in rule add-to-array action
1040				next if ! defined $rules->{$k}; # Do we update key $k?
1041				print "\t\tAdding $k with code $rules->{$k} ".
1042				"and value '$r->{value}' ...\n"
1043				if $debug =~ /h/;
1044				my $accum_is_blank = ! defined $accum->{$k} \|\|
1045				! defined $accum->{$k}{value} \|\|
1046				grep($accum->{$k}{value} eq $_, @blank);
1047				my $this_is_blank = ! defined $r->{value} \|\|
1048				grep($r->{value} eq $_, @blank);
1049				# If the rule action is an array ref, do all. The user must make sure
1050				# the actions will result in different header card names, which pretty
1051				# much means only 1 of the actions may have a 1-3 value, all others
1052				# must be 'XXX@@@@@' type rules.
1053				my @actions;
1054				if(ref $rules->{$k}) {
1055				@actions = @{$rules->{$k}};
1056				} else {
1057				@actions = ($rules->{$k});
1058				}
1059				for my $action (@actions) {
1060				if ($action eq '1') { # Keep first non-blank
1061				if($accum_is_blank) {
1062				$accum->{$k}{value} = $r->{value};
1063				}
1064				} elsif($action eq '2') { # Keep last non-blank
1065				if($accum_is_blank \|\| ! $this_is_blank) {
1066				$accum->{$k}{value} = $r->{value};
1067				}
1068				} elsif($action eq '3') { # Add non-blanks
1069				if(! $this_is_blank) {
1070				$accum->{$k}{value} += $r->{value};
1071				# Make the type always be double so we don't
1072				# get problems with precision truncating things that
1073				# look vaguely int-like. (The auto-typer in the
1074				# FITSIO module is limited in the magic it can perform.)
1075				$accum->{$k}{type} = 'F';
1076				#print "$k: new=$r->{value}, sum=".
1077				# sprintf("%22.16g",$accum->{$k}{value})."\n";
1078				}
1079				} elsif($action eq '4') { # Running avg of non-blanks
1080				if(! $this_is_blank) {
1081				++$accum->{$k}{_num};
1082				$accum->{$k}{value} +=
1083				($r->{value} - ($accum->{$k}{value}\|\|0)) /
1084				$accum->{$k}{_num};
1085				# See comment above
1086				$accum->{$k}{type} = 'F';
1087				}
1088				} elsif($action =~ /^([^\@]+)(\@+)$/) { # Append to array
1089				my ($keybase,$digs) = ($1,$2);
1090				# Construct a key name from the rule.
1091				# 'base@@@@' means make an array like base0001, base002, etc.
1092				$newname = sprintf "$keybase%0".length($digs)."d",$n;
1093				$accum->{$newname}{value} = $r->{value};
1094				$accum->{$newname}{name} = $newname;
1095				} else {
1096				confess "$0/accumhdr: Add-key code '$action' for key '$k' ".
1097				"not recognized.";
1098				}
1099				} # Actions
1100				}
1101
1102				if(! @minkeep) {
1103				@minkeep = qw/SIMPLE EXTEND XTENSION BITPIX EXTNAME
1104				NAXIS NAXIS1 NAXIS2 NAXIS3 NAXIS4
1105				BLANK GCOUNT PCOUNT COMMENT HISTORY
1106				LONGSTRN/
1107				}
1108
1109				print "NEW =\n",Dumper($accum) if $debug =~ /hh/;
1110
1111				# Establish what header cards to keep.
1112				my %keep;
1113				# Add explicitely added cards. Do it using a hash to avoid duplicates
1114				@keep{keys %$accum} = ();
1115				# Add defaults.
1116				@keep{@minkeep} = ();
1117				# Array-ize.
1118				my @keep = sort keys %keep;
1119
1120				print "Accum HDU =\n",Dumper($accum) if $debug =~ /hh/;
1121
1122				return wantarray ? ($accum,$rules,\@keep) : $accum;
1123				}
1124
1125				# Package to hold helper routines for reform_tbl
1126
1127				package FITSIO::reform;
1128
1129	1	1.0e-6	1.0e-6	my $PI = atan2(1.0,1.0) * 4.0;
1130	1	1.0e-6	1.0e-6	my $R2D = 180.0 / $PI;
1131	1	1.0e-6	1.0e-6	my $S2 = sqrt(2.0);
1132
1133				sub ecl2gal {
1134				my ($rows,$rowix,$info,@cols) = @_;
1135				require WISE::CoUtils;
1136				my @res = WISE::CoUtils::cconv('ecl',
1137				$rows->{$cols[0]\|\|'ELON'}[$rowix],
1138				$rows->{$cols[1]\|\|'ELAT'}[$rowix],
1139				'gal')
1140				or return;
1141				return @res[0,1];
1142				}
1143
1144				sub ecl2equ {
1145				my ($rows,$rowix,$info,@cols) = @_;
1146				require WISE::CoUtils;
1147				my @res = WISE::CoUtils::cconv('ecl',
1148				$rows->{$cols[0]\|\|'ELON'}[$rowix],
1149				$rows->{$cols[1]\|\|'ELAT'}[$rowix],
1150				'equ')
1151				or return;
1152				return @res[0,1];
1153				}
1154
1155				sub poserr_deg {
1156				my ($rows,$rowix,$info,@cols) = @_;
1157				my $lonerr2 = ($rows->{$cols[0]}[$rowix]*
1158				cos($rows->{$cols[2]}[$rowix]/$R2D))**2 +
1159				($rows->{$cols[1]}[$rowix]*
1160				sin($rows->{$cols[2]}[$rowix]/$R2D))**2;
1161				my $laterr2 = ($rows->{$cols[0]}[$rowix]*
1162				sin($rows->{$cols[2]}[$rowix]/$R2D))**2 +
1163				($rows->{$cols[1]}[$rowix]*
1164				cos($rows->{$cols[2]}[$rowix]/$R2D))**2;
1165				return (sqrt($lonerr2), sqrt($laterr2));
1166				}
1167
1168				sub inrownum {
1169				my ($rows,$rowix,$info,@cols) = @_;
1170				return $rowix+1;
1171				}
1172
1173				sub rownum {
1174				my ($rows,$rowix,$info,@cols) = @_;
1175				return $info->{rownum};
1176				}
1177
1178				sub bandnum {
1179				my ($rows,$rowix,$info,@cols) = @_;
1180				return WISE::Utils::bandnum($info->{band}) \|\| -1;
1181				}
1182
1183				sub bandstr {
1184				my ($rows,$rowix,$info,@cols) = @_;
1185				return WISE::Utils::bandstr($info->{band}) \|\| '';
1186				}
1187
1188				sub area_beams {
1189				my ($rows,$rowix,$info,@cols) = @_;
1190				my $beam_as = WISE::Utils::bandfwhm($info->{band}) \|\| -3600;
1191				return $rows->{$cols[0]}[$rowix]36003600/$beam_as**2
1192				}
1193
1194				sub sex2flux {
1195				my ($rows,$rowix,$info,@cols) = @_;
1196				my $pixsz = $info->{fitshdr}{CDELT1} \|\| $info->{fitshdr}{PXSCAL1}/3600;
1197				$pixsz = abs($pixsz)*3600 if $pixsz;
1198				my $conv = WISE::Utils::band_jy_pix_ta($info->{band},$pixsz)
1199				or return -999;
1200				#use WISE::Dumper;
1201				#print Dumper $info;
1202				#print "$info->{band},$pixsz,$rows->{$cols[0]}[$rowix],$conv\n";
1203				return $rows->{$cols[0]}[$rowix] * $conv * 1000; # mJy
1204				}
1205
1206	1	1.0e-5	1.0e-5	1;