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

File	/wise/base/deliv/dev/lib/perl/WISE/Wrap.pm
Statements Executed	37
Total Time	0.014314 seconds

Subroutines — ordered by exclusive time
Calls	P	F	Exclusive Time	Inclusive Time	Subroutine
0	0	0	0	0	WISE::Wrap::	BEGIN
0	0	0	0	0	WISE::Wrap::OO::	AUTOLOAD
0	0	0	0	0	WISE::Wrap::OO::	BEGIN
0	0	0	0	0	WISE::Wrap::OO::	new
0	0	0	0	0	WISE::Wrap::	_load_tbl
0	0	0	0	0	WISE::Wrap::	build_gauss_filter
0	0	0	0	0	WISE::Wrap::	clean_all
0	0	0	0	0	WISE::Wrap::	coadd_layout
0	0	0	0	0	WISE::Wrap::	gen_coadd_id
0	0	0	0	0	WISE::Wrap::	get_prex_cands
0	0	0	0	0	WISE::Wrap::	get_xsc_mchs
0	0	0	0	0	WISE::Wrap::	getfix
0	0	0	0	0	WISE::Wrap::	init_version
0	0	0	0	0	WISE::Wrap::	is_known_layout
0	0	0	0	0	WISE::Wrap::	lm_from_jhk
0	0	0	0	0	WISE::Wrap::	make_img_jpeg
0	0	0	0	0	WISE::Wrap::	make_img_mJPEG
0	0	0	0	0	WISE::Wrap::	mk_list_file
0	0	0	0	0	WISE::Wrap::	normalize_work
0	0	0	0	0	WISE::Wrap::	orbitnum_from_scanid
0	0	0	0	0	WISE::Wrap::	path_from_scanid
0	0	0	0	0	WISE::Wrap::	path_root_subst
0	0	0	0	0	WISE::Wrap::	pixint
0	0	0	0	0	WISE::Wrap::	print_jk_tbl
0	0	0	0	0	WISE::Wrap::	read_meta_tbl
0	0	0	0	0	WISE::Wrap::	run_idl
0	0	0	0	0	WISE::Wrap::	scanframe_from_path
0	0	0	0	0	WISE::Wrap::	scanframe_to_path
0	0	0	0	0	WISE::Wrap::	scangrp_from_scanid
0	0	0	0	0	WISE::Wrap::	set_bands
0	0	0	0	0	WISE::Wrap::	steps_resolve
0	0	0	0	0	WISE::Wrap::	true_percent
0	0	0	0	0	WISE::Wrap::	update_wcs_from_meta
0	0	0	0	0	WISE::Wrap::	wcs_cd_to_scalerot
0	0	0	0	0	WISE::Wrap::	wcs_cdelt_to_cd
0	0	0	0	0	WISE::Wrap::	wcs_rot
0	0	0	0	0	WISE::Wrap::	wcs_scale
0	0	0	0	0	WISE::Wrap::	wcs_scalerot
0	0	0	0	0	WISE::Wrap::	write_gaussian_filter
0	0	0	0	0	WISE::Wrap::	zexist
0	0	0	0	0	WISE::Wrap::	zlink

Line	Stmts.	Exclusive Time	Avg.	Code
1				#! /usr/bin/env perl
2
3	3	4.2e-5	1.4e-5	use strict; # spent 13µs making 1 call to strict::import
4	3	6.2e-5	2.1e-5	use warnings; # spent 34µs making 1 call to warnings::import
5
6				package WISE::Wrap;
7
8	3	3.9e-5	1.3e-5	use WISE::Env (cfglib=>'/wise/base/deliv/dev/lib/perl', import=>[qw/$cfgdir $opsdir $staticdir/]); # spent 639µs making 1 call to WISE::Env::import, max recursion depth 1
9
10				#use Carp;
11
12				# Copied from WISE::CoUtils so I don't have to include the module
13	3	6.5e-5	2.2e-5	use vars qw /$PI $R2D $S2 $sig2fwhm/; # spent 58µs making 1 call to vars::import
14	1	1.0e-6	1.0e-6	$PI = 3.1415926535897931;
15	1	1.0e-6	1.0e-6	$R2D = 180.0 / $PI;
16	1	1.0e-6	1.0e-6	$S2 = sqrt(2.0);
17	1	1.0e-6	1.0e-6	$sig2fwhm = 2.354820;
18
19	1	1.0e-6	1.0e-6	my $version = '$Id: Wrap.pm 7879 2010-05-11 23:07:01Z tim $ ';
20
21	3	2.9e-5	9.7e-6	use vars qw(@ISA @EXPORT @EXPORT_OK $VERSION); # spent 63µs making 1 call to vars::import
22
23	3	0.01376	0.00459	use Exporter::Lite; # spent 39µs making 1 call to Exporter::Lite::import
24
25	1	1.0e-6	1.0e-6	$VERSION = 1.00;
26	1	3.0e-6	3.0e-6	@ISA = ();
27
28				sub run_idl {
29				my $command = shift;
30				my $opts = shift \|\| {};
31				my $idl = $opts->{idl} \|\| 'idl';
32				my $idlerr_re = $opts->{idlerr_re};
33				my $lines = $opts->{out};
34				my $idlpath = $opts->{idlpath} \|\| '';
35				my $idldir = $opts->{idldir} \|\| '';
36				my $idlrc = $opts->{idlrc};
37				my $test = $opts->{test};
38				my $verbose = $opts->{verbose} \|\| $test;
39				my $iam = $opts->{iam} \|\| "$0";
40				my $err = "*** $iam/RUNIDL";
41				my ($pid,$in,$out);
42
43				require IPC::Open2; IPC::Open2->import();
44
45				$idlpath = join(":",@$idlpath) if ref $idlpath;
46				if($idlpath =~ s/^=//) {
47				$idlpath =~ s/^none$//i;
48				$ENV{IDL_PATH} = $idlpath;
49				} else {
50				$ENV{IDL_PATH} = (ref($cfgdir)
51				? join(",",map {"+$cfgdir/lib/idl:"} @$cfgdir)
52				: "+$cfgdir/lib/idl:").
53				"+$opsdir/lib/idl:".
54				"+$staticdir/lib/idl:".
55				"+/wise/cots/idl/idl/lib:".
56				"+/SciApps/idlastrodir/idlastro/pro:".
57				"+/home/kam/deepfields/idl:".
58				"+/home/kam/mirlin/idl".
59				($idlpath ? ":$idlpath" : "");
60				}
61				if($idlrc) {
62				$ENV{IDL_STARTUP} = ($idlrc =~ /^env/i
63				? $ENV{IDL_STARTUP}
64				: $idlrc);
65				} else {
66				$ENV{IDL_STARTUP} = "";
67				}
68				if($idldir) {
69				$ENV{IDL_DIR} = $idldir;
70				$ENV{RSI_DIR} = "$idldir/..";
71				}
72
73				print "IDL_DIR = $ENV{IDL_DIR}\n" if $verbose;
74				print "IDL_PATH = $ENV{IDL_PATH}\n" if $verbose;
75
76				print "IDL command being sent to '$idl':\n$command\n" if $verbose;
77
78				return 1 if $test;
79
80				eval { $pid = open2($in, $out, "$idl 2>&1") };
81				warn("$err: Couldn't start process '$idl'; $@\n"), return
82				if $@ =~ /^open.:/i;
83
84				print $out "$command\n";
85				close $out;
86
87				my @lines = <$in>;
88				close $in;
89
90				@$lines = @lines if ref $lines;
91
92				my @err = grep { /^\*/ \|\|
93				/^% Execution halted at/ \|\|
94				/\bERROR\b/ \|\|
95				/syntax error/i \|\|
96				/READU: End of file/ \|\|
97				($idlerr_re && /$idlerr_re/)
98				} @lines;
99				warn("$err: IDL command failed:\n$command\n", @lines), return
100				if @err;
101
102				print @lines if $verbose;
103
104				waitpid($pid,0);
105
106				warn("$err: IDL command ended with code $?:\n$command\n"), return
107				if $?;
108
109				return wantarray ? (@lines) : 1;
110				}
111
112				# Write a FITS gaussian filter file
113				sub write_gaussian_filter {
114				my $file = shift; # Output FITS file name
115				my $fwhm = shift; # Width of PSF, in pixels
116				my $n = shift; # Filter size in pixel
117				my $cdelt = shift; # Pixel width in degrees
118				my $opts = shift \|\| {};
119				my $add = $opts->{addhdr} \|\| {};
120				my @add = $add ? ( ref($add)=~/hash/i ? (%$add) : (@$add) ) : ();
121				my $err = "*** $0";
122
123				require WISE::FITSIO;
124				require PDL; PDL->import();
125
126				my $filter = build_gauss_filter($fwhm,$n,$opts);
127
128				my $fits = WISE::FITSIO->new($file,{mode=>"w",
129				hdudef=>[[
130				BITPIX =>
131				$opts->{bitpix}\|\|-64,
132				NAXES => [$n,$n],
133				CDELT2 => $cdelt,
134				FWHM => $fwhm,
135				@add,
136				]]
137				})
138				or warn("$err: Cannot open PSF output file '$file'.\n"), return;
139
140				$fits->writepix($filter->get_dataref())
141				or warn("$err: Cannot write PSF to file '$file'.\n"), return;
142
143				return wantarray ? $filter : 1;
144				}
145
146				# require PDL; PDL->import();
147				# require PDL::Math; PDL::Math->import('erf');
148				# require PDL::Image2D; PDL::Image2D->import();
149
150				# Build an $n x $n, normalized, optionally zero sum, gaussian filter
151				sub build_gauss_filter {
152				my $fwhm = shift; # Width of PSF, in pixels
153				my $n = shift \|\| int($fwhm+1)*2 + 1;
154				my $opts = shift \|\| {};
155				my $zsum = $opts->{zsum};
156				my ($x,$y,$sum);
157				my @filter;
158
159				require PDL; PDL->import();
160
161				my ($min,$max) = (int(-$n/2), ($n%2 ? int($n/2) : int($n/2)-1));
162
163				for $y ($min..$max) {
164				for $x ($min..$max) {
165				my $z = pixint(0,0,$fwhm,$x,$y);
166				push @{$filter[$y+int($n/2)]}, $z;
167				$sum += $z;
168				}
169				}
170				# Normalize
171				for my $i (0..$n-1) { for my $j (0..$n-1) {
172				$filter[$i][$j] /= $sum;
173				} }
174				if($zsum) {
175				# Make zero-sum
176				for my $i (0..$n-1) { for my $j (0..$n-1) {
177				$filter[$i][$j] -= 1/$n**2;
178				} }
179				}
180				my $pdl = ! $opts->{bitpix} \|\| $opts->{bitpix}==-64
181				? pdl(double(),\@filter)
182				: pdl(float(), \@filter);
183
184				return $pdl;
185				}
186
187				# Return the normalized flux in a unit pixel relative to a gaussian at
188				# a given center and fwhm
189				sub pixint {
190				my ($x0,$y0,$fwhm,$x,$y) = @_;
191				# $x/y0 = centroid of gaussian
192				# fwhm = FWHM of gaussian
193				# $x,$y = center of pixel to integrate
194
195				require PDL; PDL->import();
196				require PDL::Math; PDL::Math->import('erf');
197
198				my $sigma = $fwhm / $sig2fwhm;
199				my $xhi = ($x+.5-$x0)/($sigma*$S2);
200				my $xlo = ($x-.5-$x0)/($sigma*$S2);
201				my $yhi = ($y+.5-$y0)/($sigma*$S2);
202				my $ylo = ($y-.5-$y0)/($sigma*$S2);
203				my ($erxhi) = list(erf($xhi));
204				my ($erxlo) = list(erf($xlo));
205				my ($eryhi) = list(erf($yhi));
206				my ($erylo) = list(erf($ylo));
207
208				return 0.25($erxhi - $erxlo)($eryhi - $erylo);
209				}
210
211				# Stub for routine now moved to Params.pm (for compatability only)
212				sub steps_resolve {
213				require WISE::Params;
214				return WISE::Params::steps_resolve(@_);
215				}
216
217				sub read_meta_tbl {
218				my $tbl = shift;
219				my $opts = shift \|\| {};
220				my $iam = $opts->{iam} \|\| "$0/READMETA";
221				my $simple = $opts->{simple};
222				my $err = "*** $iam/META";
223				my ($data,$meta,%meta);
224				require WISE::IOUtils;
225				eval {
226				($data,$meta) = WISE::IOUtils::slurp_ipac_tbl($tbl, {lc=>1});
227				};
228				die "$err: Couldn't read meta-data table '$tbl'.\n$@\n"
229				if ! $data;
230				for (0..$meta->{nrows}-1) {
231				my $name = lc $data->{name}[$_];
232				$name =~ s/^.*:// if $opts->{short};
233				if($simple) {
234				$meta{$name}[$data->{band}[$_]] = $data->{value}[$_];
235				} else {
236				$meta{$name}[$data->{band}[$_]] = {val=>$data->{value}[$_],
237				hdr=>$data->{hdrname}[$_],
238				comment=>$data->{comment}[$_],
239				};
240				}
241				}
242
243				return wantarray ? %meta : \%meta;
244				}
245
246				# Fully breakdown a CD matrix either from a FITS filename, a FITSIO object,
247				# a header hash, or a bare CD element array.
248				sub wcs_cd_to_scalerot {
249				my $hdr = shift;
250				my $file = "<UNKNOWN>";
251				my @cd;
252				my @cdhdr;
253				die "*** $0/WCSCD2SCLROTBETA: Too many params passed.\n"
254				if @_;
255				if(ref($hdr) =~ /ARRAY/) {
256				@cdhdr = @$hdr;
257				warn("*** $0/WCSCSCLROT: CD matrix array does not have 4 elements.\n"),
258				return
259				if @cdhdr != 4;
260				@cd = (undef,[undef,$cdhdr[0],$cdhdr[1]],[undef,$cdhdr[2],$cdhdr[3]]);
261				$file = 'parameter CD array';
262				} else {
263				if(! ref $hdr) {
264				my $tmphdr = WISE::IOUtils::fits_hdr($hdr);
265				warn("*** $0/WCSCD2SCLROTBETA: Unable to get header from file '$hdr'.\n"),
266				return
267				if ! %$tmphdr;
268				$file = $hdr;
269				$hdr = $tmphdr;
270				} elsif(ref($hdr) =~ /fitsio/i) {
271				$hdr = $hdr->key(["*"],{kv=>1});
272				warn("*** $0/WCSCSCLROT: Unable to get header from file '".$hdr->{file}."'.\n"),
273				return
274				if ! $hdr;
275				$file = $hdr->{file};
276				} elsif(ref($hdr) =~ /HASH/) {
277				$file = 'parameter hash';
278				} else {
279				die "*** $0/WCSCD2SCLROTBETA: Don't know what to do with ref of type '".ref($hdr)."'.\n";
280				}
281				my @cdkeys = grep {defined $hdr->{$_}} qw/CD1_1 CD1_2 CD2_1 CD2_2/;
282				for (sort @cdkeys) {
283				next if ! defined $_;
284				my ($i,$j) = /CD(\d)_(\d)/;
285				$cd[$i][$j] = $hdr->{$_};
286				push @cdhdr, $hdr->{$_};
287				}
288				}
289				warn("*** $0/WCSCSCLROT: Unable to get 4 CD matrix elements from $file.\n"),
290				return
291				if @cdhdr != 4;
292				# From HLM brkdcd.f fortran routine
293				my @rcd;
294				my $rho = atan2(-$cd[1][2],$cd[2][2]);
295				$rcd[1][1] = cos($rho)$cd[1][1] + sin($rho)$cd[2][1];
296				$rcd[1][2] = cos($rho)$cd[1][2] + sin($rho)$cd[2][2];
297				$rcd[2][1] = -sin($rho)$cd[1][1] + cos($rho)$cd[2][1];
298				$rcd[2][2] = -sin($rho)$cd[1][2] + cos($rho)$cd[2][2];
299				my $s2 = $rcd[2][2];
300				my $beta = POSIX::asin(-$rcd[2][1]/$rcd[1][1]);
301				my $s1 = $rcd[1][1]*cos($beta);
302				my $crota2 = WISE::Utils::angnorm($rho*$R2D);
303				my $pa = WISE::Utils::angnorm(-$crota2);
304				return wantarray ? ($s1,$s2,$crota2,$beta)
305				: {cdelt1=>$s1,cdelt2=>$s2,crota2=>$crota2,
306				skew=>$beta,cd=>\@cdhdr,pa=>$pa};
307				}
308
309				sub wcs_scalerot {
310				my $hdr = shift;
311				if(ref($hdr) =~ /fitsio/i) {
312				$hdr = $hdr->key(["*"],{kv=>1});
313				warn("*** $0/WCSCSCLROT: Unable to get header from ".$hdr->{file}.".\n"),
314				return
315				if ! $hdr;
316				}
317				my $crota2 = wcs_rot($hdr);
318				my ($cdelt1,$cdelt2) = wcs_scale($hdr,$crota2);
319				return ($cdelt1,$cdelt2,$crota2);
320				}
321
322				sub wcs_rot {
323				my $hdr = shift;
324				if(ref($hdr) =~ /fitsio/i) {
325				$hdr = $hdr->key(["*"],{kv=>1});
326				warn("*** $0/WCSCROT: Unable to get header from ".$hdr->{file}.".\n"),
327				return
328				if ! $hdr;
329				}
330				my ($crota2);
331				$crota2 = $hdr->{CROTA2};
332				$crota2 = $hdr->{WCROTA2}
333				if ! defined $crota2;
334				if(! defined $crota2) {
335				(undef,undef,$crota2) = wcs_cd_to_scalerot($hdr) or return;
336				}
337				warn("$0/WCSCROT: Can't get CROTA2.\n"), return
338				if ! defined $crota2;
339				$crota2 += 360 if $crota2 < 0;
340				return $crota2;
341				}
342
343				sub wcs_scale {
344				my $hdr = shift;
345				my $crota2 = shift; # Optional
346				if(ref($hdr) =~ /fitsio/i) {
347				$hdr = $hdr->key(["*"],{kv=>1});
348				warn("*** $0/WCSSCL: Unable to get header from ".$hdr->{file}.".\n"),
349				return
350				if ! $hdr;
351				}
352				my ($cdelt1,$cdelt2);
353				($cdelt1,$cdelt2) = @{$hdr}{qw/CDELT1 CDELT2/};
354				($cdelt1,$cdelt2) = @{$hdr}{qw/WCDELT1 WCDELT2/}
355				if ! $cdelt1 \|\| ! $cdelt2;
356				if(! $cdelt1 \|\| ! $cdelt2) {
357				# Use CD matrix
358				($cdelt1,$cdelt2,$crota2) = wcs_cd_to_scalerot($hdr) or return;
359				}
360				warn("$0/WCSSCL: Can't get CDELT's.\n"), return
361				if ! $cdelt1 \|\| ! $cdelt2;
362				my $cdelt = $cdelt1 ? abs($cdelt1) : abs($cdelt2);
363				return wantarray ? ($cdelt1,$cdelt2) : $cdelt;
364				}
365
366				sub wcs_cdelt_to_cd {
367				my ($cdelt1,$cdelt2,$crota2,$skew,$opts) = @_;
368				my $err = "*** $0/WCS2CD";
369				($opts,$skew) = ($skew,undef) if ! $opts && ref $skew;
370				require POSIX;
371				if(ref($cdelt1)) { # Single header hash passed
372				my $hdr = $cdelt1;
373				($opts,$cdelt2) = ($cdelt2,undef) if ! $opts && ref $cdelt2;
374				die "$err: Unwanted args defined.\n"
375				if $cdelt2 \|\| $crota2 \|\| $skew;;
376				if(ref($hdr) =~ /fitsio/i) {
377				$hdr = $hdr->key(["*"],{kv=>1});
378				warn("*** $0/WCSCD: Unable to get header from ".
379				$hdr->{file}.".\n"),
380				return
381				if ! $hdr;
382				}
383				($cdelt1,$cdelt2) = wcs_scale($hdr);
384				($crota2,$skew) = @{$hdr}{qw/CROTA2 SKEW/};
385				$crota2 = -$hdr->{PA} if ! defined $crota2;
386				} else {
387				die "$err: Options in wrong place?\n"
388				if ref $cdelt2 \|\| ref $crota2 \|\| ref $skew;;
389				}
390				$opts \|\|= {};
391				my $more = $opts->{more};
392				my $comment = $opts->{comments};
393				$crota2 \|\|= 0;
394				my $pa = $crota2 > 0 ? 360-$crota2 : -$crota2;
395				$skew \|\|= 0;
396				my $cd11 = $cdelt1*cos($crota2/$R2D);
397				my $cd12 = -$cdelt2*sin($crota2/$R2D);
398				my $cd21 = $cdelt1*sin($crota2/$R2D);
399				my $cd22 = $cdelt2*cos($crota2/$R2D);
400				if(abs($skew) > 1e-15) {
401				$cd11 = $cd11 + -$cd12$skew$cdelt1/$cdelt2;
402				$cd21 = $cd21 + -$cd22$skew$cdelt1/$cdelt2;
403				}
404				my @new = ( CD1_1 => $cd11. ($comment?' / WCS rotation matrix element' :''),
405				CD1_2 => $cd12. ($comment?' / WCS rotation matrix element' :''),
406				CD2_1 => $cd21. ($comment?' / WCS rotation matrix element' :''),
407				CD2_2 => $cd22. ($comment?' / WCS rotation matrix element' :''),
408				($more
409				? (CDELT1 => $cdelt1. ($comment?' / [deg/pix] X-axis scale' :''),
410				CDELT2 => $cdelt2. ($comment?' / [deg/pix] Y-axis scale' :''),
411				CROTA2 => $crota2. ($comment?' / [deg] CCW rotation of RA at CRPIX1,2' :''),
412				WCROTA2 => $crota2. ($comment?' / [deg] CCW rotation of RA at CRPIX1,2' :''),
413				WCDELT1 => $cdelt1. ($comment?' / [deg/pix] X-axis scale' :''),
414				WCDELT2 => $cdelt2. ($comment?' / [deg/pix] Y-axis scale' :''),
415				SKEW => $skew. ($comment?' / [deg] Y-axis shear' :''),
416				PA => $pa. ($comment?' / [deg] Rotation of +Y EofN at CRPIX1,2' :''),
417				PXSCAL1 => ($cdelt1*3600).($comment?' / [arcsec/pixel] X-axis scale' :''),
418				PXSCAL2 => ($cdelt2*3600).($comment?' / [arcsec/pixel] Y-axis scale' :''),
419				)
420				: ()
421				)
422				);
423				return wantarray ? @new : { @new };
424				}
425
426				# Update frames with refined WCS values from a meta-file
427				sub update_wcs_from_meta {
428				my $imgfile = shift;
429				my $metafile = shift;
430				my $opts = shift \|\| {};
431				my $app = $opts->{app} \|\| "";
432				my $meta_band= $opts->{meta_band};
433				my $newlvl = $opts->{newlvl}\|\| "";
434				my $nodistort= $opts->{nodistort}; # Zero out SIP distortion header values
435				my $nrep = $opts->{nrep} \|\| ($newlvl =~ /[23]/ ? 2 : 1); # !!!!!
436				my $dup = $opts->{dupwcs}; # Don't apply cirrection, just copy file
437				my $err = "*** ".($opts->{iam}\|\|"$0/WCSMETA");
438				my $warn= "=== ".($opts->{iam}\|\|"$0/WCSMETA");
439				require WISE::UtilsLight;
440				require WISE::FITSIO;
441				# Keys/rules from meta-data table needed for updating WCS.
442				# Use array to preserve order. Derive hash below.
443				$app .= ":" if $app;
444				my @update_keys_hdrs =("${app}bfraf" => [["CRVAL1", 1]],
445				"${app}bfdcf" => [["CRVAL2", 1]],
446				"${app}bfpaf" => [["CROTA2", -1],
447				["WCROTA2", -1],
448				["PA", 1]],
449				"${app}sigraf" => [["CRDER1", 1/3600.]],
450				"${app}sigdcf" => [["CRDER2", 1/3600.]],
451				"${app}sigpaf" => [["UNCRTPA", 1/3600.]],
452				"${app}cosigrdf" => [["CSDRADEC", 1/3600]],
453				"${app}sigsxf" => [["UNCRTS1", 1]],
454				"${app}sigsyf" => [["UNCRTS2", 1]],
455				(! $newlvl \|\| $newlvl =~ /1/
456				? ("${app}bfsfxf" => [["CDELT1", 1./$nrep],
457				["WCDELT1", 1./$nrep],
458				["PXSCAL1", 3600./$nrep]],
459				"${app}bfsfyf" => [["CDELT2", 1./$nrep],
460				["WCDELT2", 1./$nrep],
461				["PXSCAL2", 3600./$nrep]])
462				# Don't want to alter scales for coadds
463				: ()
464				)
465				);
466				my %update_keys_hdrs = @update_keys_hdrs; # When order not important
467				my @update_keys = WISE::UtilsLight::aryhash_keys(@update_keys_hdrs);
468				# WCS values needed for computing the CD array
469				my @cdkeys = ("${app}bfsfxf", "${app}bfsfyf", "${app}bfpaf");
470
471				my $realfile = zexist($imgfile);
472
473				die "$err: Image file '$imgfile' not found.\n"
474				if ! $realfile;
475
476				# Get WCS correction meta data
477				my %meta;
478				if($metafile && ! ref $metafile) {
479				print "Reading meta-table '$metafile' ...\n" if $opts->{verbose};
480				die "$err: Meta-data table '$metafile' not found.\n"
481				if ! -e $metafile;
482				%meta = read_meta_tbl($metafile,{short=>$app?0:1, simple=>1});
483				} elsif($metafile && keys %$metafile) {
484				%meta = %$metafile;
485				} else {
486				# Missing PRex data is OK if we're dup'ing
487				if(! $dup) {
488				die "$err: No PRex meta data found.\n";
489				}
490				}
491				my $parts = WISE::UtilsLight::pathdecomp($realfile);
492				$meta_band \|\|= WISE::UtilsLight::bandnum($parts->{band});
493				my $fr = $realfile;
494				my $newfr = defined $newlvl
495				? WISE::UtilsLight::pathcomp($realfile,{vsn=>$newlvl})
496				: $fr;
497				if(! WISE::UtilsLight::samefile($fr,$newfr)) {
498				print "Upgrading WCS of '$fr' to '$newfr' ...\n" if $opts->{verbose};
499				if(! $opts->{overwrite}) {
500				require File::Copy;
501				print "Copying '$fr' to '$newfr' ...\n" if $opts->{verbose};
502				File::Copy::copy($fr,$newfr)
503				or die "$err: Failed to copy '$fr' to '$newfr'; $!.\n";
504				} else {
505				print "Renaming '$fr' to '$newfr' ...\n" if $opts->{verbose};
506				rename($fr,$newfr)
507				or die "$err: Failed to rename '$fr' to '$newfr'; $!.\n";
508				}
509				} else {
510				print "Upgrading WCS of '$newfr' in place ...\n" if $opts->{verbose};
511				}
512
513				my $fits = WISE::FITSIO->new($newfr,{mode=>'rw',verbose=>$opts->{verbose}})
514				or die "$err: Unable to open frame for update; file='$newfr'.\n";
515
516				if($dup) {
517				warn "$warn: Copying WCS data directly from source image unchanged.\n";
518				return $fits;
519				}
520
521				# Assemble new WCS keys
522				my (@new);
523				for my $key (@update_keys) {
524				die "$err: No meta-data found for '$key', band $meta_band in ".
525				"'$metafile'.\n"
526				if ! defined $meta{$key}[$meta_band];
527				for my $ref (@{$update_keys_hdrs{$key}}) {
528				my ($hdrkey,$scale) = @$ref;
529				push @new, ($hdrkey => $scale*$meta{$key}[$meta_band]);
530				}
531				}
532				# Compute new CD_ key values (and others)
533				# (Use of aryhash_keys() allows me to treat an array of k,v pairs
534				# like a hash, but preserve order.)
535				if($newlvl =~ /1/) { # Don't want a CD matrix for coadds
536				my ($cdelt1, $cdelt2, $pa) = map {$meta{$_}[$meta_band]} @cdkeys;
537				if($cdelt1 == 0 \|\| $cdelt2 == 0) {
538				# !!! Later this should die
539				warn "$err: No valid meta-data found for band $meta_band in ".
540				"'$metafile'. Use original data.\n";
541				return 1;
542				}
543				my $skew = $fits->key("SKEW") \|\| 0;
544				push @new, WISE::Wrap::wcs_cdelt_to_cd($cdelt1,$cdelt2,-$pa,$skew,
545				{
546				more=>1, # WCDELT/PXSCAL/PA too
547				comments=>1, # Add HDU comments
548				}
549				);
550				}
551				# Save old as comments
552				my @old = $fits->key([WISE::UtilsLight::aryhash_keys(@new)], {kv=>1});
553				# Filter out WCS values not present (many will not be in level 2/3)
554				{
555				my %old = @old;
556				@old = map { ($_=>$old{$_}) }
557				grep {defined $old{$_}}
558				WISE::UtilsLight::aryhash_keys(@old);
559				my %new = @new;
560				@new = map { ($_=>$new{$_}) }
561				grep {defined $old{$_}}
562				WISE::UtilsLight::aryhash_keys(@new);
563				}
564				# Update with new
565				print "Updating WCS values: '@new' ...\n" if $opts->{verbose};
566				$fits->key([\@new]);
567				if($nodistort) {
568				print "!!! Zeroing out distortion ...\n" if $opts->{verbose};
569				my @distort = map { ($_=>0) } grep {/^[AB]P?(_\d+_\d+\|_DMAX)$/}
570				WISE::UtilsLight::aryhash_keys($fits->keyhash(['*']));
571				if(@distort) { $fits->key([\@distort]); }
572				}
573				print "Saving old WCS values: '@old' ...\n" if $opts->{verbose};
574				$fits->key([["COMMENT" => "OLD ".shift(@old)." = ".shift(@old)]])
575				while @old;
576				# (Leaving any distortion model unchanged.)
577
578				return $fits;
579				}
580
581
582				sub make_img_jpeg {
583				my $imgfile = shift;
584				my $opts = shift \|\| {};
585				my $band = $opts->{band};
586				require WISE::UtilsLight;
587				my $jpegfile = $opts->{jpegfile}
588				\|\| WISE::UtilsLight::pathcomp($imgfile,{form=>'jpg'});
589				my $regfiles = $opts->{regfiles} \|\| [];
590
591				my $parts = WISE::UtilsLight::pathdecomp($imgfile);
592				$band \|\|= $parts->{band};
593
594				print "Initiating ds9 session ...\n" if $opts->{verbose};
595				my $ds9 = WISE::DS9->new({verbose=>$opts->{verbose},
596				blank=>1,tile=>1,min=>1});
597				$ds9->sz(900);
598				my $fr = $imgfile;
599				print "Displaying frame '$fr' ...\n" if $opts->{verbose};
600				$ds9->disp($fr,{new=>1,zoom=>$band&&$band==4?1:0.5,mode=>"99.5"});
601				-e $_ and $ds9->reg($_) for @{$regfiles};
602				$ds9->grid();
603				$ds9->fr(1);
604				$ds9->align();
605				$ds9->mch();
606				print "Writing jpeg '$jpegfile' ...\n" if $opts->{verbose};
607				$ds9->jpeg($jpegfile);
608
609				return 1;
610				}
611
612
613				sub make_img_mJPEG {
614				my $opts = shift \|\| {};
615				my $iam = $opts->{iam} \|\| "$0";
616				my $err = "*** $iam/mJPEG";
617				require WISE::Spawn;
618
619				my $stretch = $opts->{logpower} \|\| 'log';
620				my $min = $opts->{min} \|\| '0%';
621				my $max = $opts->{max} \|\| '98.5%';
622
623				my @args = ('mJPEG');
624				if($opts->{gray}) {
625				my ($mymin,$mymax) = ($min,$max);
626				if($opts->{mypct} && "$min$max"=~/%/) {
627				($mymin,$mymax) = true_percent($opts->{gray},$min,$max);
628				print "Generated new gray min,max = $mymin,$mymax from $min,$max\n"
629				if $opts->{verbose};
630				}
631				push @args, ( '-gray', $opts->{gray}, $mymin, $mymax, $stretch);
632				} elsif( $opts->{red} && $opts->{green} && $opts->{blue} ) {
633				my ($mymin_r,$mymax_r,$mymin_g,$mymax_g,$mymin_b,$mymax_b) = ($min,$max,$min,$max,$min,$max);
634				if($opts->{mypct} && "$min$max"=~/%/) {
635				($mymin_r,$mymax_r) = true_percent($opts->{red},$min,$max);
636				print "Generated new red min,max = $mymin_r,$mymax_r from $min,$max\n"
637				if $opts->{verbose};
638				($mymin_g,$mymax_g) = true_percent($opts->{green},$min,$max);
639				print "Generated new green min,max = $mymin_g,$mymax_g from $min,$max\n"
640				if $opts->{verbose};
641				($mymin_b,$mymax_b) = true_percent($opts->{blue},$min,$max);
642				print "Generated new green min,max = $mymin_b,$mymax_b from $min,$max\n"
643				if $opts->{verbose};
644				}
645				push @args, ( '-red', $opts->{red}, $min, $max, $stretch,
646				'-green', $opts->{green}, $min, $max, $stretch,
647				'-blue', $opts->{blue}, $min, $max, $stretch );
648				} else {
649				warn "$err: incorrect inputs. Specify either \"gray\" or \"red\", \"green\", and \"blue\". \n";
650				return;
651				}
652				push @args, ( '-out', $opts->{jpegfile}, '-norotate');
653
654				my $rc = WISE::Spawn::spawncmd(\@args,{verbose=>$opts->{verbose},
655				test=>$opts->{test},
656				noexit=>1});
657				if($rc) {
658				warn "$err: mJPEG return value $rc"; return;
659				}
660
661				return 1;
662				}
663
664				sub true_percent {
665				my $img = shift;
666				my @pcts = @_;
667				my $err = "*** $0/true_percent";
668				require WISE::FITSIO::PDL;
669				my $i = WISE::FITSIO::PDL->new($img);
670				if(! $i) {
671				warn "$err: Unable to open '$img'.\n";
672				return;
673				}
674				my $pdl = $i->getpdl();
675				my @x;
676				for my $p (@pcts) {
677				if($p =~ s/%//) {
678				my $ptile = $pdl->pct($p/100.);
679				if($ptile =~ /nan/) { $ptile = $p; } # Punt
680				push @x, $ptile;
681				} else {
682				push @x, $p;
683				}
684				}
685				return @x;
686				}
687
688
689				sub path_from_scanid {
690				my $frame_id = shift;
691				my $opts = shift \|\| {};
692				my $top = $opts->{pipedir} \|\| "/wise/fops/scans";
693				my ($scan_id,$frnum,$src) = $frame_id =~
694				m%
695				(?:^\|[\s/])
696				(\d{5}[a-z]) # Scan ID
697				(\d{3}) # Frame num
698				(?:[-_\s/.](\d{6})\|$)
699				%x;
700				return if ! $scan_id;
701				my $scangrp = scangrp_from_scanid($scan_id);
702				my $dir = "$top/$scangrp/$scan_id/fr/$frnum";
703				return wantarray ? ($dir,$scan_id.$frnum,$scan_id,$frnum,$src) : $dir;
704				}
705
706				sub scanframe_from_path {
707				my $path = shift;
708				return if ! $path;
709				my ($scanid,$framenum,$scangrp,$root,$base,$file);
710				my ($file_scanid,$file_scangrp,$file_framenum);
711				my $warn = "=== $0/SCANIDFROMPATH";
712				require WISE::UtilsLight;
713
714				if($path =~ m\|/\|) {
715				# Only resolve if a path is attached
716				my $tmp = eval {
717				WISE::UtilsLight::normalizepath($path,{resolve=>1,notrail=>1});
718				};
719				$tmp = WISE::UtilsLight::normalizepath($path) if ! $tmp; # No abs_path()
720				$path = $tmp \|\| $path;
721				}
722				($root,$scangrp,$scanid,$framenum,$file) =
723				$path =~ m%^
724				(.*) # root
725				/+
726				(\d[\da-z]) # scan group
727				/+
728				(\d{5}[a-z]) # scan ID
729				(?: /+fr/+
730				(\d\d\d) # frame number
731				)? # ... is optional
732				(?:
733				/+([^/]*) # trailing file name
734				)? # ... is optional
735				$
736				%x;
737				if($scanid) {
738				# Got info from dir, compare to file
739				($file_scanid,$file_framenum) =
740				$file =~ m%^
741				(\d\d\d\d\d[a-z]) # scan ID
742				(\d\d\d) # frame number
743				(?:[-.]\|$)
744				%x
745				if $file && $framenum;
746				if($file_scanid) {
747				$file_scangrp = scangrp_from_scanid($file_scanid);
748				# File also matched; compare results
749				warn "$warn: File base ($file_scanid$file_framenum) ".
750				"and directory ($scanid$framenum) scan/frame components ".
751				"disagree for path '$path'.\n"
752				if $file_scanid.$file_framenum ne $scanid.$framenum;
753				}
754				} else {
755				# Dir not matched; try to get info from file name
756				($scanid,$framenum) =
757				$path =~ m%^(?:.*/+)? # optional directiry
758				(\d\d\d\d\d[a-z]) # scan ID
759				(\d\d\d) # frame number
760				(?:[-.]\|$) # termination
761				%x;
762				$scangrp = scangrp_from_scanid($scanid) if $scanid;
763				}
764
765				if($scanid && $framenum) {
766				$base = $scanid.$framenum;
767				}
768				return if ! $scanid; # Couldn't match anything useful
769				my $orbitnum = orbitnum_from_scanid($scanid);
770				return (wantarray
771				? ($scanid,$framenum,$scangrp,$root,$base,$orbitnum)
772				: {scanid=>$scanid,framenum=>$framenum,scangrp=>$scangrp,
773				root=>$root,base=>$base,orbitnum=>$orbitnum}
774				);
775				}
776
777				sub scanframe_to_path {
778				my ($scanid,$framenum,$scangrp,$root,$base) = @_;
779				require WISE::UtilsLight;
780				if(ref $scanid) {
781				my %tmp = %$scanid;
782				($scanid,$framenum,$scangrp,$root,$base) =
783				@tmp{qw/scanid framenum scangrp root base/};
784				}
785				return if ! $scanid;
786				$scangrp //= scangrp_from_scanid($scanid);
787				my $path = sprintf("%s/%s".
788				(defined $framenum ? "/fr/%03d" : "").
789				"/%s",
790				$scangrp,$scanid,
791				(defined $framenum ? ($framenum) : ()),
792				$base\|\|"");
793				$path = "$root/$path" if $root;
794				return WISE::UtilsLight::normalizepath($path,{notrail=>1});
795				}
796
797				sub scangrp_from_scanid {
798				return substr($_[0]\|\|"00000x",4,2);
799				}
800
801				sub orbitnum_from_scanid {
802				return substr($_[0]\|\|"00000x",0,5) / 2;
803				}
804
805				sub get_prex_cands {
806				my ($ra0,$dec0) = (shift,shift);
807				my ($opts,$image);
808				$image = $ra0, $opts = $dec0 \|\| {}, $ra0 = $dec0 = undef
809				if ! defined $dec0 \|\| ref($dec0);
810				$opts = shift \|\| {}
811				if ! $opts;
812				my $refdb = $opts->{refdb} \|\| "/wise/fops/ref/astrom/2mass_ref.db";
813				my $reftbl = $opts->{reftbl} \|\| "ref_2mass";
814				my $refinfo = $opts->{refinfo} \|\| $reftbl."_info";
815				my $r = $opts->{r} \|\| 0.4; # FOV inner circle
816				my $noradtest = $opts->{noradtest};
817				my $maxrows = $opts->{maxrows};
818				my $behead = $opts->{behead};
819				my $outtbl = $opts->{out};
820				my $ds9 = $opts->{ds9};
821				my $j_k_limit = $opts->{j_k_limit} \|\| 0.5;
822				my $l_limit = $opts->{l_limit} \|\| 99;
823				my $harsh = $opts->{harsh};
824				my $retrymin = $opts->{retrymin} \|\| 10;
825				my $ds9color = $opts->{ds9color} \|\| 'cyan';;
826				my $err = "*** ".($opts->{iam}\|\|"$0/GETPREX");
827				my $warn = "=== ".($opts->{iam}\|\|"$0/GETPREX");
828				my ($racol,$deccol) = qw(ra dec);
829				my %type_cvt = (
830				"double" => 'r',
831				"real" => 'r',
832				"float" => 'r',
833				"int" => 'i',
834				"char" => 'c',
835				);
836				my %col_cvt = (
837				decl => 'dec'
838				);
839
840				eval "use DBI; 1;"
841				or die "$err: Can't load DBI.\n$@";
842				eval "use WISE::CHealPix; 1;"
843				or die "$err: Can't load WISE::CHealPix.\n$@";
844				eval "use WISE::Time; 1;"
845				or die "$err: Can't load WISE::Time.\n$@";
846				eval "use WISE::CoUtils; 1;"
847				or die "$err: Can't load WISE::CoUtils.\n$@";
848
849				if($image) {
850				# Get ra,dec0 from image header
851				eval "use WISE::FITSIO; 1;"
852				or die "$err: Can't load WISE::CoUtils.\n$@";
853				my $f;
854				if(! ref $image) {
855				$f = WISE::FITSIO->new($image) or die;
856				} else {
857				$f = $image;
858				$image = $f->file();
859				}
860				print "Getting center positions from '$image' ...\n"
861				if $opts->{verbose};
862				my $h = $f->key(["*"],{kv=>1});
863				($ra0,$dec0) = @{$h}{qw/CRVAL1 CRVAL2/};
864				die "$err: CRVAL1,2 keywords not present in '$image'.\n"
865				if ! defined $ra0 \|\| ! defined $dec0;
866				}
867
868				print "Connecting to $refdb ...\n" if $opts->{verbose};
869				my $d=DBI->connect("DBI:SQLite:dbname=$refdb","","",{RaiseError=>1});
870				my $nside;
871				NSIDE: {
872				my $s = $d->prepare("select nside from $refinfo");
873				$s->execute;
874				($nside) = $s->fetchrow_array();
875				die "$err: Cannot extract NSIDE from '$refinfo'.\n"
876				if ! $nside;
877				}
878				print "Read NSIDE=$nside from '$refinfo' table.\n" if $opts->{verbose};
879				my $pix0 = WISE::CHealPix::ang2pix_ring($nside,$ra0,$dec0);
880				my $newr = $r + WISE::CHealPix::c2cdist_deg($nside)/sqrt(2);
881				my @hpbins = WISE::CHealPix::pix_withinr_ring($nside,$pix0,$newr);
882				print "Found ".@hpbins." bins near $ra0,$dec0 (pix=$pix0, r=$newr).\n"
883				if $opts->{verbose};
884				my $harsh_query = $harsh
885				? "AND j_m - k_m < $j_k_limit \n"
886				: "";
887				my ($sel,$nrows,$ngood);
888				my (@cols,@sqltypes,@types,@lens,@blanks,@sortmag);
889				my %rows;
890				RERUN_QUERY: {
891				my $sel = "SELECT * FROM $reftbl WHERE hpbin in (".
892				join(",",@hpbins).") \n".
893				$harsh_query.
894				"ORDER BY k_m";
895				print "Filtering with following selection:\n$sel\n"
896				if $opts->{verbose};
897				$nrows = 0;
898				$ngood = 0;
899				(@cols,@sqltypes,@types,@lens,@sortmag,%rows) = ();
900				QUERY: {
901				my $s = $d->prepare($sel);
902				$s->execute;
903				@cols = map { $col_cvt{$_}\|\|$_ } @{$s->{NAME_lc}};
904				@sqltypes = @{$s->{TYPE}};
905				@blanks = ('null') x @cols;
906				for my $t (@sqltypes) {
907				push @types, (map {$type_cvt{$_}}
908				grep {$t=~/$_/}
909				keys %type_cvt)
910				[0];
911				}
912				print "COLS = @cols\n" if $opts->{debug};
913				print "SQLTYPES = @sqltypes\n" if $opts->{debug};
914				print "TYPES = @types\n" if $opts->{debug};
915				die "$err: Some columns did not match a known type pattern.\n"
916				if grep {! $_} @types;
917				my %colnums = map { ($cols[$_]=>$_) } 0..$#cols;
918				my ($rano,$decno) = @colnums{$racol, $deccol};
919				my ($jno,$hno,$kno) = @colnums{qw/j_m h_m k_m/};
920				die "$err: Cannot find columns $racol,$deccol in table $reftbl.\n"
921				if ! defined $racol \|\| ! defined $deccol;
922				FETCH: while(my $row=$s->fetchrow_arrayref()) {
923				++$nrows;
924				next if $noradtest \|\|
925				WISE::CoUtils::angdist($ra0,$dec0,
926				$row->[$rano],$row->[$decno])
927				> $r;
928				if($harsh_query) {
929				my ($l,$m) = lm_from_jhk($row->[$jno],
930				$row->[$hno],
931				$row->[$kno]);
932				next if ! defined $l;
933				next if $l > $l_limit;
934				push @sortmag, $l;
935				} else {
936				push @sortmag, $row->[$kno];
937				}
938				++$ngood;
939				for my $c (0..$#cols) {
940				push @{$rows{$cols[$c]}}, $row->[$c];
941				my $l = length($row->[$c]);
942				$lens[$c] = ! defined $lens[$c] \|\| $l>$lens[$c]-2
943				? $l+2
944				: $lens[$c];
945				}
946				if($opts->{debug} && $opts->{debug}>1) {
947				print $cols[$_],'=',$rows{$cols[$_]}[-1],',' for 0..$#cols;
948				print "\n";
949				}
950				} # FETCH
951				} # QUERY
952				if($harsh_query && $ngood < $retrymin) {
953				$harsh_query = ''; # Must be an empty string (not '0' etc.)
954				print "Re-running query less harshly since number usable ($ngood) ".
955				"is < $retrymin ...\n"
956				if $opts->{verbose};
957				redo RERUN_QUERY;
958				} else {
959				$d->disconnect;
960				$d = undef;
961				}
962				print "Query returned $nrows rows, of which $ngood will be used.\n"
963				if $opts->{verbose};
964				} # RERUN_QUERY
965				warn "$warn: Query returned no sources < $r degrees from $ra0,$dec0.\n"
966				if $ngood == 0;
967				# Add sortmag column
968				push @cols, 'sort_m';
969				push @types, 'r';
970				push @lens, '7';
971				push @blanks, 'null';
972				push @sqltypes, '';
973				$rows{sort_m} = \@sortmag;
974				if($harsh_query) {
975				# Sort on l instead of k
976				my @sortix = sort {$sortmag[$a]<=>$sortmag[$b]} 0..$#sortmag;
977				@{$rows{$_}} = @{$rows{$_}}[@sortix] for keys %rows;
978				}
979				if($behead && $maxrows > 0 && $ngood > $behead) {
980				# Remove top rows
981				@{$rows{$_}} = @{$rows{$_}}[$behead..$ngood-1] for keys %rows;
982				$ngood -= $behead;
983				print "Reducing row count to $ngood after removing top $behead sources.\n"
984				if $opts->{verbose};
985				}
986				if($maxrows && $maxrows > 0 && $ngood > $maxrows) {
987				# Reduce to desired row count
988				@{$rows{$_}} = @{$rows{$_}}[0..$maxrows-1] for keys %rows;
989				$ngood = $maxrows;
990				print "Reducing row count to $ngood brightest.\n"
991				if $opts->{verbose};
992				}
993				$rows{_nrows} = $ngood;
994				$rows{_cols} = \@cols;
995				$rows{_sqltypes} = \@sqltypes;
996				$rows{_types} = \@types;
997				$rows{_lens} = \@lens;
998				if($outtbl && $ngood > 0) {
999				eval "use WISE::IOUtils; 1;"
1000				or die "$err: Can't load WISE::IOUtils.\n$@";
1001				my $datime = WISE::Time::Time_str(time(),{form=>1});
1002				print "Writing table '$outtbl' ...\n" if $opts->{verbose};
1003				WISE::IOUtils::make_ipac_tbl($outtbl\|\|"-",
1004				{
1005				names =>\@cols,
1006				types =>\@types,
1007				lens =>\@lens,
1008				blanks=>\@blanks,
1009				},
1010				[
1011				{name=>RowsRetreived => value=>$ngood },
1012				{name=>RowsRetrieved => value=>$ngood },
1013				{name=>NRows => value=>$ngood },
1014				{name=>DateTime => value=>$datime},
1015				{name=>DBName => value=>$refdb },
1016				{name=>DBTable => value=>$reftbl},
1017				{name=>NSide => value=>$nside },
1018				{name=>RA0 => value=>$ra0 },
1019				{name=>Dec0 => value=>$dec0 },
1020				{name=>SearchRadius => value=>$r },
1021				],
1022				\%rows,
1023				{
1024				lc=>1,
1025				})
1026				or die "$err: Could not write table '$outtbl'.\n";
1027				}
1028				if($ds9 && $ngood > 0) {
1029				print "Writing ds9 regions to '$ds9' ...\n" if $opts->{verbose};
1030				WISE::IOUtils::write_ds9_regions($ds9,\%rows,
1031				{
1032				coord => 'j2000',
1033				loncol => 'ra',
1034				latcol => 'dec',
1035				color => $ds9color,
1036				symbol => 'cross',
1037				a0 => 5/3600,
1038				})
1039				or warn "$warn: Failed to write regions to '$ds9'.\n";
1040				}
1041
1042				return wantarray ? %rows : \%rows;
1043				}
1044
1045				# Try to generate WISE-like w1 and w2 fluzes if possible
1046
1047	1	1.0e-6	1.0e-6	my ($minjk,$maxjk,%jk_tbl);
1048	1	3.0e-6	3.0e-6	my @jk_cols = qw/V-K J-K J-H H-K K-L V-12 K-M sigma SpT/;
1049	1	1.0e-6	1.0e-6	my $djk = 0.01;
1050
1051				sub lm_from_jhk {
1052				my ($j,$h,$k) = @_;
1053				my $jk = sprintf("%6.2f",$j-$k);
1054				_load_tbl() if !%jk_tbl;
1055				return if $jk < $minjk \|\| $jk > $maxjk;
1056				return if ! defined $jk_tbl{$jk}{'K-L'} \|\| ! defined $jk_tbl{$jk}{'K-M'};
1057				return (
1058				$k-$jk_tbl{$jk}{'K-L'}, # L
1059				$k-$jk_tbl{$jk}{'K-M'} # M
1060				);
1061				}
1062
1063				sub print_jk_tbl {
1064				_load_tbl() if !%jk_tbl;
1065				printf +("%6s " x 8)."%6s\n",@jk_cols;
1066				for my $jk (sort {$a<=>$b} keys %jk_tbl) {
1067				printf +("%6.2f " x 8)."%6s\n",@{$jk_tbl{$jk}}{@jk_cols};
1068				}
1069				}
1070
1071				sub _load_tbl {
1072				my @rawtbl;
1073				# Read rows
1074				while(<DATA>) {
1075				next if /^#/ \|\| /^\s*$/; # Skip blank and comment liness
1076				push @rawtbl,[split " "];
1077				my %row;
1078				@row{@jk_cols} = @{$rawtbl[-1]};
1079				$minjk = $row{'J-K'} if ! defined $minjk;
1080				$maxjk = $row{'J-K'};
1081				}
1082				# Interpolate on 0.01 j-k values
1083				for my $rownum (0..$#rawtbl-1) {
1084				my (%row0,%row1);
1085				@row0{@jk_cols} = @{$rawtbl[$rownum]};
1086				@row1{@jk_cols} = @{$rawtbl[$rownum+1]};
1087				my $maxdjk = $row1{'J-K'} - $row0{'J-K'};
1088				my $maxdkl = $row1{'K-L'} - $row0{'K-L'};
1089				my $maxdkm = $row1{'K-M'} - $row0{'K-M'};
1090				my $ndjk = int( $maxdjk/$djk ) + 1;
1091				for my $jk_ix (0..$ndjk) {
1092				my $deljk= $jk_ix*$djk;
1093				my $jk = sprintf("%6.2f",$row0{'J-K'} + $deljk);
1094				my $kl = $row0{'K-L'} + ($deljk/$maxdjk)*$maxdkl;
1095				my $km = $row0{'K-M'} + ($deljk/$maxdjk)*$maxdkm;
1096				# ... interpolate other columns as desired ...
1097				my %row = %row0;
1098				$row{'J-K'} = $jk;
1099				$row{'K-L'} = $kl;
1100				$row{'K-M'} = $km;
1101				$jk_tbl{$jk} = \%row;
1102				}
1103				}
1104				return 1;
1105				}
1106
1107				sub get_xsc_mchs {
1108				my ($ra0,$dec0) = (shift,shift);
1109				my ($opts,$image);
1110				$image = $ra0, $opts = $dec0 \|\| {}, $ra0 = $dec0 = undef
1111				if ! defined $dec0 \|\| ref($dec0);
1112				$opts = shift \|\| {}
1113				if ! $opts;
1114				my $refdb = $opts->{refdb} \|\| "/wise/fops/ref/astrom/2mass_xref.db";
1115				my $reftbl = $opts->{reftbl} \|\| "xref_2mass";
1116				my $refinfo = $opts->{refinfo} \|\| $reftbl."_info";
1117				my $r = $opts->{r} \|\| 0.4; # FOV inner circle
1118				my $noradtest = $opts->{noradtest};
1119				my $maxrows = $opts->{maxrows};
1120				my $outtbl = $opts->{out};
1121				my $ds9 = $opts->{ds9};
1122				my $harsh = $opts->{harsh};
1123				my $retrymin = $opts->{retrymin} \|\| 10;
1124				my $ds9color = $opts->{ds9color} \|\| 'magenta';;
1125				my $err = "*** ".($opts->{iam}\|\|"$0/GETXSC");
1126				my $warn = "=== ".($opts->{iam}\|\|"$0/GETXSC");
1127				my ($racol,$deccol) = qw(ra dec);
1128				my %type_cvt = (
1129				"double" => 'r',
1130				"real" => 'r',
1131				"float" => 'r',
1132				"int" => 'i',
1133				"char" => 'c',
1134				);
1135				my %col_cvt = (
1136				decl => 'dec'
1137				);
1138
1139				eval "use DBI; 1;"
1140				or die "$err: Can't load DBI.\n$@";
1141				eval "use WISE::CHealPix; 1;"
1142				or die "$err: Can't load WISE::CHealPix.\n$@";
1143				eval "use WISE::Time; 1;"
1144				or die "$err: Can't load WISE::Time.\n$@";
1145				eval "use WISE::CoUtils; 1;"
1146				or die "$err: Can't load WISE::CoUtils.\n$@";
1147
1148				if($image) {
1149				# Get ra,dec0 from image header
1150				eval "use WISE::FITSIO; 1;"
1151				or die "$err: Can't load WISE::CoUtils.\n$@";
1152				my $f;
1153				if(! ref $image) {
1154				$f = WISE::FITSIO->new($image) or die;
1155				} else {
1156				$f = $image;
1157				$image = $f->file();
1158				}
1159				print "Getting center positions from '$image' ...\n"
1160				if $opts->{verbose};
1161				my $h = $f->key(["*"],{kv=>1});
1162				($ra0,$dec0) = @{$h}{qw/CRVAL1 CRVAL2/};
1163				die "$err: CRVAL1,2 keywords not present in '$image'.\n"
1164				if ! defined $ra0 \|\| ! defined $dec0;
1165				}
1166
1167				print "Connecting to $refdb ...\n" if $opts->{verbose};
1168				my $d=DBI->connect("DBI:SQLite:dbname=$refdb","","",{RaiseError=>1});
1169				my $nside;
1170				NSIDE: {
1171				my $s = $d->prepare("select nside from $refinfo");
1172				$s->execute;
1173				($nside) = $s->fetchrow_array();
1174				die "$err: Cannot extract NSIDE from '$refinfo'.\n"
1175				if ! $nside;
1176				}
1177				print "Read NSIDE=$nside from '$refinfo' table.\n" if $opts->{verbose};
1178				my $pix0 = WISE::CHealPix::ang2pix_ring($nside,$ra0,$dec0);
1179				my $newr = $r + WISE::CHealPix::c2cdist_deg($nside)/sqrt(2);
1180				my @hpbins = WISE::CHealPix::pix_withinr_ring($nside,$pix0,$newr);
1181				print "Found ".@hpbins." bins near $ra0,$dec0 (pix=$pix0, r=$newr).\n"
1182				if $opts->{verbose};
1183				my $harsh_query = $harsh
1184				? ""
1185				: "";
1186				my ($sel,$nrows,$ngood);
1187				my (@cols,@sqltypes,@types,@lens,@sortmag);
1188				my %rows;
1189				RERUN_QUERY: {
1190				my $sel = "SELECT * FROM $reftbl WHERE hpbin in (".
1191				join(",",@hpbins).") \n".
1192				$harsh_query."\n".
1193				"ORDER BY k_m_ext";;
1194				print "Filtering with following selection:\n'$sel'\n"
1195				if $opts->{verbose};
1196				$nrows = 0;
1197				$ngood = 0;
1198				(@cols,@sqltypes,@types,@lens,@sortmag,%rows) = ();
1199				QUERY: {
1200				my $s = $d->prepare($sel);
1201				$s->execute;
1202				@cols = map { $col_cvt{$_}\|\|$_ } @{$s->{NAME_lc}};
1203				@sqltypes = @{$s->{TYPE}};
1204				for my $t (@sqltypes) {
1205				push @types, (map {$type_cvt{$_}}
1206				grep {$t=~/$_/}
1207				keys %type_cvt)
1208				[0];
1209				}
1210				print "COLS = @cols\n" if $opts->{debug};
1211				print "SQLTYPES = @sqltypes\n" if $opts->{debug};
1212				print "TYPES = @types\n" if $opts->{debug};
1213				die "$err: Some columns did not match a known type pattern.\n"
1214				if grep {! $_} @types;
1215				my %colnums = map { ($cols[$_]=>$_) } 0..$#cols;
1216				my ($rano,$decno) = @colnums{$racol, $deccol};
1217				die "$err: Cannot find columns $racol,$deccol in table $reftbl.\n"
1218				if ! defined $racol \|\| ! defined $deccol;
1219				my ($jno,$hno,$kno) = @colnums{qw/j_m_ext h_m_ext k_m_ext/};
1220				FETCH: while(my $row=$s->fetchrow_arrayref()) {
1221				++$nrows;
1222				$_ = undef for grep { ! defined $_ \|\| /^\s*$/ } @$row;
1223				next if ! $noradtest &&
1224				(! defined $row->[$rano] \|\| ! defined $row->[$decno]);
1225				next if $noradtest \|\|
1226				WISE::CoUtils::angdist($ra0,$dec0,
1227				$row->[$rano],$row->[$decno])
1228				> $r;
1229				if($harsh_query) {
1230				push @sortmag, $row->[$kno];
1231				} else {
1232				push @sortmag, $row->[$kno];
1233				}
1234				++$ngood;
1235				for my $c (0..$#cols) {
1236				push @{$rows{$cols[$c]}}, $row->[$c];
1237				my $l = length($row->[$c]//'');
1238				$row->[$c] = undef if $l == 0;
1239				$lens[$c] = ! defined $lens[$c] \|\| $l>$lens[$c]-2
1240				? $l+2
1241				: $lens[$c];
1242				}
1243				if($opts->{debug} && $opts->{debug}>1) {
1244				print $cols[$_],'=',$rows{$cols[$_]}[-1],',' for 0..$#cols;
1245				print "\n";
1246				}
1247				last FETCH if $maxrows && $ngood >= $maxrows;
1248				} # FETCH
1249				} # QUERY
1250				$d->disconnect;
1251				print "Query returned $nrows rows, of which $ngood will be used.\n"
1252				if $opts->{verbose};
1253				if($harsh_query && $ngood < $retrymin) {
1254				$harsh_query = ''; # Must be an empty string (not '0' etc.)
1255				print "Re-running query less harshly since number usable ".
1256				"is < $retrymin ...\n"
1257				if $opts->{verbose};
1258				redo RERUN_QUERY;
1259				}
1260				} # RERUN_QUERY
1261				warn "$warn: Query returned no sources < $r degrees from $ra0,$dec0.\n"
1262				if $ngood == 0;
1263				# Add sortmag column
1264				push @cols, 'sort_m';
1265				push @types, 'r';
1266				push @lens, '7';
1267				push @sqltypes, '';
1268				$rows{sort_m} = \@sortmag;
1269				$rows{_nrows} = $ngood;
1270				$rows{_cols} = \@cols;
1271				$rows{_sqltypes} = \@sqltypes;
1272				$rows{_types} = \@types;
1273				$rows{_lens} = \@lens;
1274				if($outtbl && $ngood > 0) {
1275				eval "use WISE::IOUtils; 1;"
1276				or die "$err: Can't load WISE::IOUtils.\n$@";
1277				my $datime = WISE::Time::Time_str(time(),{form=>1});
1278				print "Writing table '$outtbl' ...\n" if $opts->{verbose};
1279				WISE::IOUtils::make_ipac_tbl($outtbl\|\|"-",
1280				{
1281				names =>\@cols,
1282				types =>\@types,
1283				lens =>\@lens,
1284				blanks=>[('null') x @cols],
1285				},
1286				[
1287				{name=>RowsRetreived => value=>$ngood },
1288				{name=>RowsRetrieved => value=>$ngood },
1289				{name=>NRows => value=>$ngood },
1290				{name=>DateTime => value=>$datime},
1291				{name=>DBName => value=>$refdb },
1292				{name=>DBTable => value=>$reftbl},
1293				{name=>NSide => value=>$nside },
1294				{name=>RA0 => value=>$ra0 },
1295				{name=>Dec0 => value=>$dec0 },
1296				{name=>SearchRadius => value=>$r },
1297				],
1298				\%rows,
1299				{
1300				lc=>1,
1301				})
1302				or die "$err: Could not write table '$outtbl'.\n";
1303				}
1304				if($ds9 && $ngood > 0) {
1305				print "Writing ds9 regions to '$ds9' ...\n" if $opts->{verbose};
1306				WISE::IOUtils::write_ds9_regions($ds9,\%rows,
1307				{
1308				coord => 'j2000',
1309				loncol => 'ra',
1310				latcol => 'dec',
1311				color => $ds9color,
1312				symbol => 'box',
1313				a0 => 5/3600,
1314				})
1315				or warn "$warn: Failed to write regions to '$ds9'.\n";
1316				}
1317
1318				return wantarray ? %rows : \%rows;
1319				}
1320
1321				# Thin wrapper for running getfix application.
1322				# Arguably this is inside out; There should be a true API that the
1323				# getfix app uses, rather than vice-versa.
1324				sub getfix {
1325				my $params = shift;
1326				my $opts = shift \|\| {};
1327				my $file = $opts->{file} // "-";
1328				my $live = $opts->{livefix};
1329				my $tblio_opts = $opts->{tblio_opts};
1330				my $defs = $opts->{defs} // $WISE::Env::pars;
1331				my $exec = 'getfix';
1332				my $err = "*** $0/getfix";
1333				my $warn = "=== $0/getfix";
1334				require WISE::Spawn;
1335				require WISE::IPACTbl;
1336				my $spawn = {
1337				take => [qw/data_root/],
1338				add => {
1339				%$params,
1340				out_file => $file,
1341				verbose => 0,
1342				(! $live
1343				? (frame_index_file =>
1344				"%frame_index_dir%/%frame_index_base%-fix-ro.db")
1345				: ()),
1346				},
1347				noexit => 1,
1348				test => 0,
1349				verbose=> $opts->{verbose},
1350				err => \ my $errout,
1351				out => \ my $out,
1352				tee => 0,
1353				};
1354				my $rc;
1355				if(! $defs) { $rc = WISE::Spawn::spawncmd($exec,$spawn); }
1356				else { $rc = $defs->spawn($exec,$spawn); }
1357				if($rc) {
1358				warn @$errout if $errout && @$errout;
1359				warn "$err: Getfix failed with code $rc.\n";
1360				return;
1361				}
1362				my $tbl;
1363				if($file eq "-") {
1364				if(! $out) {
1365				warn "$warn: No table output.\n";
1366				return;
1367				}
1368				$file = \$out;
1369				}
1370				return 1 if $opts->{noread};
1371				$tbl = WISE::IPACTbl->new($file,$tblio_opts)
1372				or die "$err: Unable to read getfix output from '$file'.\n";
1373				my $rows = $tbl->data();
1374				if(! $rows) {
1375				warn "$warn: Unable to instantiate table rows from '$file'.\n";
1376				return;
1377				}
1378				return wantarray ? ($rows,$tbl) : $rows;
1379				}
1380
1381				# Take a list of files or directories and make a file of list
1382				# for a certain type of file out of it
1383
1384				sub mk_list_file {
1385				my $file = shift;
1386				my $parts= shift \|\| {};
1387				my $list = shift;
1388				my $opts = shift \|\| {};
1389				my $dirs = $opts->{dirs};
1390				my $missingok = $opts->{missingok};
1391				my $err = "*** $0/MKLIST";
1392				my $warn= "=== $0/MKLIST";
1393				require WISE::UtilsLight;
1394				my @out;
1395
1396				my $fh;
1397				if($file && ! ref $file) {
1398				open $fh, ">", $file
1399				or die "$err/list: Unable to open frame list file '$file'; $!.\n";
1400				} elsif(ref $file) {
1401				$fh = $file;
1402				$file = "<Unknown>";
1403				}
1404				if(! ref $list) {
1405				# $list if an extant list file from which we wish to retrieve
1406				# the file/dir list
1407				open my $listfh, "<", $list
1408				or die "$err: Unable to open frame list file '$list'; $!.\n";
1409				my @frames = (
1410				grep {! /^$/ && ! /^#/ }
1411				map {chomp; s/^\s//; s/\s$//; $_;}
1412				<$listfh>
1413				);
1414				$list = \@frames;
1415				}
1416				for my $f (@$list) {
1417				my $base;
1418				my $isdir = 0;
1419				if(-d $f \|\| $dirs) {
1420				my ($scanid,$framenum) = WISE::Wrap::scanframe_from_path($f);
1421				die "$err: Unable to determine base name from frame ".
1422				"directory '$f'.\n"
1423				if ! $scanid \|\| ! $framenum;
1424				$base = $scanid.$framenum;
1425				$isdir = 1;
1426				}
1427				my $new = WISE::UtilsLight::pathcomp($f,
1428				{
1429				dir => $parts->{dir},
1430				base => $base,
1431				band => $parts->{band},
1432				type => $parts->{type},
1433				vsn => $parts->{vsn},
1434				form => 'fits'
1435				},
1436				{
1437				isdir => $isdir,
1438				}
1439				);
1440
1441				if(! -e $new && ($opts->{noz} \|\| ! -e "$new.gz")) {
1442				die "$err: Requested file '$new' (from '$f', isdir=$isdir) ".
1443				"does not exist.\n"
1444				if ! $missingok;
1445				warn "$warn: Requested file '$new' (from '$f', isdir=$isdir) ".
1446				"does not exist.\n";
1447				next;
1448				}
1449
1450				if($fh) {
1451				print $fh "$new\n"
1452				or die "$err: Unable to write to list file '$file'; $!.\n";
1453				}
1454
1455				push @out, $new;
1456				}
1457
1458				if($fh) {
1459				close $fh
1460				or die "$err: Unable to write to (close) list file '$file'; $!.\n";
1461				}
1462
1463				return wantarray ? @out : \@out;
1464				}
1465
1466				# Clean a specific set of standard directories according to a standard rule set
1467
1468				sub clean_all {
1469				my $dir = shift;
1470				my $parts= shift \|\| {};
1471				my $opts = shift \|\| {};
1472				my $indir = $opts->{in_dir} // $dir; # For clean_dir safety
1473				my $outdir = $opts->{out_dir} // $dir; # "
1474				my $workdir= $opts->{work_dir} // "$dir/work";
1475				my $locworkdir= $opts->{local_work_dir} // "$dir/work";
1476				my $qadir = $opts->{qa_dir} // "$dir/qa";
1477				my $caldir = $opts->{cal_dir} // "$dir/cal";
1478				my $err = "*** $0/CLEANALL";
1479				my $warn = "*** $0/CLEANALL";
1480				require WISE::IOUtils;
1481				require Cwd;
1482
1483				return 1 if $parts->{none};
1484
1485				print +($opts->{recursing}?"Recursively c":"C")."leaning '$dir', ".
1486				"parts=".join(",",grep {$parts->{$_}} keys %$parts)." ".
1487				"from CWD=".Cwd::getcwd()." ...\n"
1488				if $opts->{verbose};
1489
1490				print "Looking at versions ...\n" if $opts->{debug};
1491
1492				if($parts->{all} \|\| $parts->{vers}) {
1493				for my $ver (grep {-d $_ && m\|$dir/v\d+$\|} glob("$dir/*")) {
1494				print "Cleaning version directory '$ver' ...\n"
1495				if $opts->{verbose};
1496				clean_all($ver,{%$parts, vers=>0}, {top_dir_too=>1,
1497				recursing=>1,
1498				verbose=>$opts->{verbose},
1499				});
1500				}
1501				}
1502
1503				print "Looking at work dir '$locworkdir' ...\n" if $opts->{debug};
1504
1505				if($parts->{all} \|\| $parts->{work}) {
1506				# Clean the outdir 'work' dir, which may be a symlink
1507				if(-d $locworkdir \|\| -l $locworkdir) {
1508				print "Cleaning work directory '$locworkdir' ".
1509				(-l $locworkdir ? "(->'".readlink($locworkdir)."') " : "").
1510				"...\n"
1511				if $opts->{verbose};
1512				# Extra warning to check if the automounter is sometimes failing to
1513				# recognize access attempts to a dir through a symlink
1514				warn "$warn: Work dir link '$locworkdir' to '".
1515				readlink($locworkdir)."' ".
1516				"fails -d check.\n"
1517				if ! -d $locworkdir && -l $locworkdir;
1518				# This will often be a symlink to node-local storage
1519				WISE::IOUtils::clean_dir($locworkdir,{dir_too=>1,
1520				parents_too=>1,
1521				recurse=>1,
1522				trash=>1,
1523				trashdir=>'default',
1524				in_dir =>$indir,
1525				out_dir=>$outdir,
1526				verbose=>$opts->{verbose},
1527				})
1528				or return;
1529				}
1530				# Also pre-clean the new work dir which may pre-exist.
1531				# Might be different from above.
1532				if(-d $workdir) {
1533				print "Pre-clean new work dir '$workdir' (?) ".
1534				(-l $workdir ? "(->'".readlink($workdir)."') " : "").
1535				"...\n"
1536				if $opts->{verbose};
1537				WISE::IOUtils::clean_dir($workdir,{in_dir =>$indir,
1538				out_dir=>$outdir,
1539				verbose=>$opts->{verbose},
1540				})
1541				or return;
1542				}
1543				}
1544
1545				print "Looking at qa dir '$qadir' ...\n" if $opts->{debug};
1546
1547				if(($parts->{all} \|\| $parts->{qa}) && -d $qadir) {
1548				print "Cleaning qa directory '$qadir' ...\n" if $opts->{verbose};
1549				WISE::IOUtils::clean_dir($qadir, {dir_too=>1,
1550				trash=>1,
1551				trashdir=>'default',
1552				in_dir =>$indir,
1553				out_dir=>$outdir,
1554				robust=>1,
1555				verbose=>$opts->{verbose},
1556				})
1557				or return;
1558				}
1559
1560				print "Looking at cal dir '$caldir' ...\n" if $opts->{debug};
1561
1562				if($parts->{cal} && -d $caldir) {
1563				print "Cleaning cal directory '$caldir' ...\n" if $opts->{verbose};
1564				WISE::IOUtils::clean_dir($caldir, {dir_too=>1,
1565				trash=>1,
1566				trashdir=>'default',
1567				in_dir =>$indir,
1568				out_dir=>$outdir,
1569				robust=>1,
1570				verbose=>$opts->{verbose},
1571				})
1572				or return;
1573				}
1574
1575				print "Looking at top-level dir '$dir' ...\n" if $opts->{debug};
1576
1577				if($parts->{all} \|\| $parts->{top}) {
1578				print "Cleaning top level directory '$dir' ...\n" if $opts->{verbose};
1579				WISE::IOUtils::clean_dir("$dir",
1580				{
1581				dir_too=>$opts->{top_dir_too},
1582				in_dir =>$indir,
1583				out_dir=>$outdir,
1584				robust=>1,
1585				not_matching=>
1586				['.*-w\d-int-0\.fits(\.gz)?',
1587				'.-meta-ingest.\.tbl',
1588				'.*\.log',
1589				'.*\.save_log',
1590				'.*\.condor_log',
1591				'fr',],
1592				verbose=>$opts->{verbose},
1593				})
1594				or return;
1595				}
1596				return 1;
1597				}
1598
1599				# Produce standard work directory from parameters and update parameters
1600
1601				sub path_root_subst {
1602				my $path = shift;
1603				my $match= shift;
1604				my $repl = shift;
1605				my $opts = shift \|\| {};
1606				my $resolve = $opts->{resolve};
1607				$path = WISE::UtilsLight::normalizepath($path,
1608				{resolve=>1,
1609				cwd=>$opts->{cwd}})
1610				if $resolve;
1611				s\|//+\|/\|g for ($path,$match,$repl);
1612				s\|/+$\|/\| for ($path);
1613				s\|/+$\|\| for ($match,$repl);
1614				s\|^/wise-ops/\d\d/wise/\|/wise/\| for ($path,$match,$repl);
1615				if($match !~ m\|/compute\|) {
1616				$path =~ s\|^/compute/wcnode\d+/wise/\|/wise/\|;
1617				}
1618				$path =~ s\|^$match/\|$repl/\|;
1619				return $path;
1620				}
1621
1622				sub normalize_work {
1623				my $pars = shift;
1624				my $steps = shift;
1625				my %pvals = $pars->vals();
1626				my $err = "*** $0/NORMWORK";
1627				my $warn= "=== $0/NORMWORK";
1628				require WISE::UtilsLight;
1629				require Cwd;
1630
1631				my $absout = eval { Cwd::abs_path($pvals{out_dir}); };
1632				die "$err: Unable to resolve out directory '$pvals{out_dir}'; $!.\n"
1633				if ! $absout;
1634				my $origwork = $pvals{work_dir};
1635				my $abswork = eval { Cwd::abs_path($pvals{work_dir}); };
1636				if(! $abswork && $pvals{work_dir} =~ m\|^/+local/\| &&
1637				$pvals{out_dir} !~ m\|^/+local/\|) {
1638				# Make a heuristical correction to the work directory reference.
1639				# It sometimes happens that a rerun is done from wsspipe that
1640				# want the work directory to be on /local, as though it was the
1641				# first time through. This wrong for rerun, so this is an attempt
1642				# to carry on as though the user had specified the global rather
1643				# than local work dir in wsspipe.
1644				warn "$warn: Unable to resolve local work directory ".
1645				"'$pvals{work_dir}'.\n".
1646				"$warn: Replacing local work directory with ".
1647				"'$pvals{out_dir}/work'.\n";
1648				$pvals{work_dir} = "$pvals{out_dir}/work";
1649				$abswork = eval { Cwd::abs_path($pvals{work_dir}); };
1650				}
1651
1652				die "$err: Unable to resolve work directory '$pvals{work_dir}'; $!.\n"
1653				if ! $abswork;
1654
1655				$absout =~ s\|//+\|/\|g;
1656				$abswork=~ s\|//+\|/\|g;
1657				my ($localwork,$uniformwork) = ($abswork,$abswork);
1658				# !!! This is klugey and strikes me as fragile
1659				# (No penalty for the below doing nothing)
1660				$uniformwork = path_root_subst($uniformwork,
1661				$pvals{local_root},
1662				$pvals{uniform_local_root});
1663				# Only try to use the /local work dir if it already exists
1664				## I don't think we ever really want to do this. We only want to
1665				## Use a /local/... path if it was actually specified
1666				## explicitly as the work dir, as wsspipe does when running
1667				## the frames.
1668				##$localwork = path_root_subst($localwork,
1669				## $pvals{uniform_local_root},
1670
1671				# Update the param data structure so the fully resloved work dir
1672				# will inherit properly, always using the localwork representation
1673				$pars->set("work_dir",$localwork);
1674
1675				# The 'frame_dir' is often set to be the work_dir
1676				my $absframe = eval { Cwd::abs_path($pvals{frame_dir}); }
1677				if $pvals{frame_dir};
1678				if($pvals{frame_dir} && ! $absframe &&
1679				$pvals{frame_dir} =~ m\|^/+local/./work/$\|) {
1680				$absframe = $localwork;
1681				$pars->set("frame_dir",$absframe);
1682				}
1683
1684				# Make sure a 'work' dir or link exists from out_dir as long as
1685				# . The work dir is really in /local..., and
1686				# . The work dir is not the same as the out dir.
1687				my $outwork = "$absout/work";
1688				#print "$uniformwork\n$outwork\n$absout\n$localwork\n";
1689				if(($steps && $steps->{setup}) &&
1690				! WISE::UtilsLight::samefile($uniformwork, $outwork) &&
1691				! WISE::UtilsLight::samefile($absout, $localwork) ) {
1692				print "Adding 'work' symlink to out dir ...\n";
1693				# Will die if '$outwork' already exists.
1694				WISE::UtilsLight::mysymlink($uniformwork, $outwork,
1695				{verbose=>$pvals{verbose},
1696				serious=>1, # Die on error
1697				replace=>0, # Do not overwrite; fail
1698				});
1699				if(! WISE::UtilsLight::samefile($absout, $localwork) ) {
1700				print "Adding '.from' symlink to work dir ...\n";
1701				WISE::UtilsLight::mysymlink($absout,"$uniformwork/.from",
1702				{verbose=>$pvals{verbose},
1703				serious=>1, # Die on error
1704				replace=>1, # Overwrite existing
1705				});
1706				}
1707				}
1708
1709				return ($abswork,$localwork,$uniformwork,$absframe,$outwork);
1710				}
1711
1712				# Do standard version setup from comandline params
1713
1714				sub init_version {
1715				my $pars = shift;
1716				my %pvals = $pars->vals();
1717				my $err = "*** $0/INITVSN";
1718				my $run_version;
1719
1720				if($pvals{new_version}) {
1721				$run_version = WISE::Utils::new_version($pvals{in_dir})
1722				or die "$err: Couldn't make a new version in '$pvals{in_dir}'.\n";
1723				print "Made new version '$run_version' in '$pvals{in_dir}'.\n";
1724				$pars->set("in_dir","$pvals{in_dir}/v$run_version");
1725				} else {
1726				$run_version = WISE::Utils::current_version($pvals{in_dir});
1727				if(defined $run_version) {
1728				print "Using version '$run_version' in '$pvals{in_dir}'.\n";
1729				$pars->set("in_dir","$pvals{in_dir}/v$run_version");
1730				}
1731				}
1732				if(defined $run_version) {
1733				$pars->expand_tags(); # Re-expand tags to propagate new directory
1734				%pvals = $pars->vals();
1735				}
1736
1737				# Make subsidiary directories
1738
1739				# QA
1740				WISE::Utils::mymkpath($pvals{qa_dir}, { verbose=>$pvals{verbose} })
1741				or die "$err: Unable to mkdir '$pvals{qa_dir}'; $!.\n";
1742
1743				# Cal
1744				WISE::Utils::mymkpath($pvals{local_cal_dir},{ verbose=>$pvals{verbose} })
1745				or die "$err: Unable to mkdir '$pvals{local_cal_dir}'; $!.\n";
1746
1747				# Work
1748				WISE::Utils::mymkpath($pvals{work_dir}, { verbose=>$pvals{verbose} })
1749				or die "$err: Unable to mkdir '$pvals{work_dir}'; $!.\n";
1750
1751				return $run_version;
1752				}
1753
1754				sub set_bands {
1755				my $pars = shift;
1756				my $bands_hash = shift;
1757				my $pvals_hash = shift;
1758				return if ! $pars && ! $pvals_hash;
1759				my %pvals= $pvals_hash ? %$pvals_hash : $pars->vals();
1760				my @bands= @{$pvals{bands} \|\| []};
1761				my (@gotbands,@gotfiles);
1762				for my $band (1..4) {
1763				my $fr = WISE::Utils::pathcomp({dir =>$pvals{out_dir},
1764				base=>$pvals{in_base},
1765				type=>'int',
1766				vsn =>'0',
1767				form=>"fits",
1768				band=>$band,
1769				});
1770				my $l0file = WISE::Wrap::zexist($fr);
1771				if($l0file) {
1772				push @gotbands,$band;
1773				push @gotfiles,$l0file;
1774				}
1775				}
1776				my %gotbands = map { ($_=>1) } @gotbands;
1777				my (%bands,$given);
1778				if(! @bands) {
1779				# Choose bands dynamically based on presence of L0 data
1780				@bands = sort {$a<=>$b} @gotbands;
1781				if($pars && $pars->has('bands')) {
1782				# Update param object if -bands is a param.
1783				$pars->set(bands => \@gotbands);
1784				$pars->expand_tags(); # Re-expand tags to propagate bands
1785				%pvals = $pars->vals(); # Incorporate changes from above sub calls
1786				}
1787				} else {
1788				@bands = sort {$a<=>$b} @bands;
1789				$given = 1;
1790				}
1791				$pvals{bands}= [@bands];
1792				%bands = map { ($_=>1) } @bands;
1793				%$bands_hash = %bands if $bands_hash;
1794				%$pvals_hash = %pvals if $pvals_hash;
1795				return wantarray ? @bands : {bands => \@bands,
1796				bands_hash => \%bands,
1797				l0bands => \@gotbands,
1798				l0bands_hash => \%gotbands,
1799				l0files => \@gotfiles,
1800				pvals => \%pvals,
1801				given => $given};
1802				}
1803
1804				sub zexist {
1805				my @ret;
1806				for my $f (@_) {
1807				(my $ftmp = $f) =~ s/\.gz$//;
1808				my ($fx) = grep {-e} ($ftmp,"$ftmp.gz");
1809				push @ret, $fx if $fx;
1810				}
1811				return wantarray ? @ret : $ret[0];
1812				}
1813				sub zlink {
1814				my @ret;
1815				for my $f (@_) {
1816				my ($fx) = grep {-l} ($f,"$f.gz");
1817				push @ret, $fx if $fx;
1818				}
1819				return wantarray ? @ret : $ret[0];
1820				}
1821
1822	1	1.0e-6	1.0e-6	my %layout_tbls;
1823	1	1.0e-6	1.0e-6	my %known_layouts = (a => 1);
1824
1825				sub is_known_layout {
1826				my $layout_id = shift // '';
1827				return $known_layouts{$layout_id};
1828				}
1829
1830				sub coadd_layout {
1831				my $layout = shift;
1832				my $opts = shift \|\| {};
1833				my $verbose = $opts->{verbose};
1834				my $coadd_pos = $opts->{coadd_pos};
1835				my $layout_pos = $opts->{layout_pos};
1836				my $index = $opts->{index};
1837				my $layout_file = $opts->{layout_file} \|\| "/wise/fops/ref/coadd/layouts";
1838				my $err = "*** $0/coadd_layout";
1839				require WISE::IPACTbl;
1840				require WISE::CoUtils;
1841				my %layout;
1842
1843				if($layout && $known_layouts{$layout}) {
1844				if($layout_file && -d $layout_file) {
1845				$layout_file = "$layout_file/$layout/coadd_layout-$layout.tbl";
1846				}
1847				die "$err: Coadd layout file '$layout_file' not present.\n"
1848				if ! -e $layout_file;
1849				$layout{file} = $layout_file;
1850				my $layout_tbl;
1851				my $layout_data;
1852				if(! $layout_tbls{$layout}) {
1853				$layout_tbl = WISE::IPACTbl->new($layout_file,{fast=>1}) or die;
1854				$layout_data = $layout_tbl->data({hashrow=>1});
1855				$layout_tbls{$layout} = $layout_data;
1856				} else {
1857				$layout_data = $layout_tbls{$layout};
1858				}
1859				my $layout_row;
1860				if(defined $index) {
1861				$layout_row = $layout_data->[$index];
1862				} elsif(defined $layout_pos &&
1863				defined $layout_pos->{ra} && defined $layout_pos->{dec}) {
1864				# Lookup by nearest RA,Dec
1865				my $ra0 = $layout_pos->{ra};
1866				my $dec0 = $layout_pos->{dec};
1867				my ($minr);
1868				my $i = 0;
1869				for my $row (@$layout_data) {
1870				my $r = WISE::CoUtils::angdist($row->{ra},$row->{dec},
1871				$ra0,$dec0);
1872				if(! defined $minr \|\| $r < $minr) {
1873				$minr = $r;
1874				$layout_row = $row;
1875				$index = $i;
1876				}
1877				++$i;
1878				}
1879				} else {
1880				die "$err: Insufficient info to find coadd layout.\n";
1881				}
1882				die "$err: No coadd layout found in '$layout_file'.\n"
1883				if ! $layout_row;
1884				# Use requested layout geometry, but allow Cl override
1885				$layout{ra} = $coadd_pos->{ra} // $layout_row->{ra};
1886				$layout{dec} = $coadd_pos->{dec}//$layout_row->{dec};
1887				$layout{pa} = $coadd_pos->{pa} // $layout_row->{pa}
1888				// $layout_tbl->key('pa');
1889				$layout{wx} = $coadd_pos->{sz} // $layout_tbl->key('xwidth')
1890				// 1.564444;
1891				$layout{wy} = $coadd_pos->{sz} // $layout_tbl->key('ywidth')
1892				// 1.564444;
1893				# Not really a standard layout if any part overridden
1894				$layout = 'z' if $coadd_pos && (defined $coadd_pos->{ra} \|\|
1895				defined $coadd_pos->{dec} \|\|
1896				(defined $coadd_pos->{sz} &&
1897				$coadd_pos->{sz} != $layout{wx}) \|\|
1898				defined $coadd_pos->{pa});
1899				} else {
1900				$layout{ra} = $coadd_pos->{ra};
1901				$layout{dec} = $coadd_pos->{dec};
1902				$layout{pa} = $coadd_pos->{pa};
1903				$layout{wx} = $coadd_pos->{sz} // 1.564444;
1904				$layout{wy} = $coadd_pos->{sz} // 1.564444;
1905				$layout \|\|= 'z';
1906				$index = -1;
1907				}
1908
1909				$layout{pa} //= 0;
1910				$layout{id} = $layout;
1911				$layout{index} = $index;
1912
1913				{
1914				my @nope = grep {! defined $layout{$_}} qw/ra dec pa wx wy/;
1915				die "$err: incomplete layout. Missing @nope\n"
1916				if @nope;
1917				}
1918
1919				return wantarray ? %layout : \%layout;
1920				}
1921
1922				sub gen_coadd_id {
1923				my ($pos,$layout_id,$tail,$opts) = @_;
1924				$opts \|\|= {};
1925				my $index = $opts->{index};
1926				my %layout = %{ $opts->{layout} \|\| {} };
1927				my $verbose = $opts->{verbose};
1928				my $use_layout = $opts->{use_layout_id};
1929				my $err = "*** $0/coadd_layout";
1930
1931				require POSIX; POSIX->import("floor");
1932
1933				if(! %layout) {
1934				%layout = WISE::Wrap::coadd_layout(
1935				$layout_id,
1936				{
1937				verbose => $verbose,
1938				($layout_id && $known_layouts{$layout_id}
1939				? (layout_pos => $pos)
1940				: (coadd_pos => $pos)
1941				),
1942				index => $index,
1943				layout_file => $opts->{layout_file},
1944				}
1945				);
1946				}
1947
1948				$layout_id= $layout{id};
1949				$index = $layout{index};
1950
1951				if(! defined $layout_id \|\| ! defined $layout{ra}) {
1952				warn "$err: Insufficient info to construct coadd ID.\n";
1953				return;
1954				}
1955
1956				#print "Constructing coadd ID ...\n" if $verbose;
1957
1958				my $rapart = sprintf("%04d", floor($layout{ra}*10));
1959				my $decpart= sprintf("%+04d",floor($layout{dec}*10));
1960				$decpart =~ s/\+/p/;
1961				$decpart =~ s/-/m/;
1962
1963				if($tail && lc($tail) ne 'none') {
1964				$tail =~ s/^[- _]+//;
1965				$tail =~ s/[- _]+$//;
1966				$tail =~ s/ +/_/g;
1967				$tail =~ s/_+/_/g;
1968				$tail =~ s/-+/_/g;
1969				} else {
1970				if(! defined $tail \|\| ! length $tail) {
1971				warn "$err: Coadd ID tail not supplied or explicitely defeated.\n";
1972				return;
1973				}
1974				$tail = undef;
1975				}
1976
1977				if($use_layout && $layout_id) {
1978				$tail = $layout_id.($tail?"_$tail":"");
1979				}
1980
1981				my $cid = $rapart.$decpart;
1982
1983				$cid .= "_$tail" if $tail;
1984
1985				if(length($tail//'') > ($opts->{tail_max}//8)) {
1986				warn "$err: Coadd tail '$tail' exceeds allowed max. length.\n";
1987				return;
1988				}
1989
1990				if(length($cid//'') > ($opts->{cid_max}//17)) {
1991				warn "$err: Coadd ID '$cid' exceeds allowed max. length.\n";
1992				return;
1993				}
1994
1995				$layout{cid} = $cid;
1996				$layout{rapart} = $rapart;
1997				$layout{ragrp} = substr($rapart,0,2);
1998				$layout{decpart} = $decpart;
1999				$layout{tail} = $tail;
2000				$layout{subdir} = "$layout{ragrp}/$layout{rapart}/".
2001				"$layout{decpart}/".
2002				($tail ? "$tail" : "");
2003
2004				return wantarray ? %layout : \%layout;
2005				}
2006
2007				#######################
2008
2009				package WISE::Wrap::OO;
2010
2011	3	0.00019	6.2e-5	use vars qw/$AUTOLOAD/; # spent 60µs making 1 call to vars::import
2012
2013				sub new {
2014				my $this = shift;
2015				my $class = ref($this) \|\| $this;
2016				return bless {},$class;
2017				}
2018
2019				# Auto-gen methods
2020				sub AUTOLOAD {
2021				my $self = shift;
2022				my $err = "*** $0/".__PACKAGE__."/AUTOLOAD";
2023				my $this = ref($self)
2024				or die "$err: '$self' is not an object.\n";
2025				return if ($AUTOLOAD =~ /::DESTROY$/);
2026				# Separate package qualifier from desired sub name
2027				my ($pkg,$sub) = $AUTOLOAD =~ m/(.:)(.)/;
2028				#print "'$AUTOLOAD'/'$pkg'/'$sub'/@_\n";
2029				# Strip off OO trailer
2030				$pkg =~ s/::OO//;
2031				# Get sub ref (do not store in namespace)
2032				my $subref;
2033				{
2034	3	0.00010	3.4e-5	no strict qw{refs}; # spent 30µs making 1 call to strict::unimport
2035				$subref = eval "\\&$pkg$sub";
2036				die "$err: Can't eval '$subref'.\n$@" if $@;
2037				}
2038				# Call
2039				goto &$subref;
2040				}
2041
2042				package WISE::Wrap;
2043
2044	1	1.4e-5	1.4e-5	1;
2045
2046				__DATA__
2047
2048				# Intrinsic colors for main sequence/giants;
2049				# adapted from Koornneef; Bessell; Jarrett
2050				# V-K J-K J-H H-K K-L V-12 K-M sigma Sp. T.
2051				0.75 -0.10 -0.15 -0.05 -0.03 -0.85 -0.08 1.0 OB
2052				0.58 -0.11 -0.07 -0.04 -0.02 -0.58 -0.07 0.4 B3
2053				0.44 -0.08 -0.06 -0.02 -0.02 -0.47 -0.05 0.4 B5
2054				0.15 -0.02 -0.01 -0.01 0.00 -0.04 0.00 0.4 A0
2055				0.08 0.03 0.02 0.00 0.00 0.05 0.00 0.4 A0
2056				0.40 0.10 0.08 0.02 0.02 0.35 -0.01 0.4 A3
2057				1.02 0.26 0.20 0.06 0.02 1.09 -0.03 0.4 F5
2058				1.17 0.31 0.24 0.07 0.04 1.49 -0.03 0.5 G0
2059				1.51 0.38 0.30 0.08 0.04 1.78 -0.04 0.5 G5
2060				2.09 0.52 0.43 0.09 0.05 2.30 -0.04 0.5 K0
2061				2.69 0.65 0.56 0.09 0.07 2.38 -0.04 0.5 K3
2062				3.20 0.73 0.62 0.10 0.10 2.95 -0.04 0.5 K5
2063				3.63 0.75 0.64 0.11 0.13 3.44 -0.05 0.5 M0
2064				4.07 0.80 0.65 0.15 0.15 3.73 -0.05 0.5 M2
2065				4.57 0.84 0.64 0.21 0.16 3.88 -0.05 0.5 M3
2066				4.98 0.88 0.63 0.25 0.16 3.92 -0.05 0.5 M4
2067				5.38 0.90 0.62 0.28 0.17 3.99 -0.05 0.6 M4
2068				5.85 0.93 0.62 0.31 0.17 4.32 -0.06 0.8 M5
2069				6.37 0.97 0.63 0.34 0.18 4.68 -0.06 1.0 M5
2070				6.97 1.01 0.65 0.36 0.18 5.20 -0.07 1.0 M6
2071				7.57 1.04 0.68 0.39 0.19 5.77 -0.07 1.0 M6
2072				8.16 1.16 0.72 0.44 0.19 6.00 -0.07 1.0 M7
2073				8.80 1.43 0.83 0.60 0.50 7.00 -0.10 0.5 L01
2074				9.35 1.71 1.02 0.70 0.63 7.00 -0.10 0.5 L23
2075				9.90 1.91 1.18 0.70 0.75 7.00 -0.10 0.5 L56
2076				10.10 1.96 1.24 0.80 0.90 7.00 -0.10 0.5 L8
2077				# giants
2078				# 2.08 0.56 0.45 0.11 0.06 2.01 -0.04 1.0 G2
2079				# 2.10 0.57 0.45 0.11 0.07 2.13 -0.04 1.0 G5
2080				# 2.16 0.59 0.47 0.12 0.07 2.27 -0.05 0.8 G8
2081				# 2.49 0.68 0.54 0.14 0.08 2.51 -0.07 0.7 K1
2082				# 2.92 0.80 0.64 0.16 0.08 3.05 -0.08 0.7 K3
2083				# 3.67 0.96 0.76 0.20 0.10 3.70 -0.09 0.6 K5
2084				# 3.74 0.97 0.76 0.21 0.12 4.08 -0.09 0.6 M0
2085				# 3.90 1.01 0.80 0.21 0.12 4.27 -0.09 0.6 M1
2086				# 4.63 1.10 0.85 0.25 0.14 4.56 -0.09 0.6 M3
2087				# 5.34 1.16 0.89 0.27 0.14 4.93 -0.09 0.6 M4
2088				# 6.20 1.24 0.93 0.31 0.15 6.00 -0.08 0.5 M5
2089				# 7.20 1.50 1.09 0.42 0.17 6.1 -0.08 0.5 M7
2090				# 9.00 3.00 1.50 1.50 1.00 10.0 1.00 1.0 AGB