| File | /wise/base/deliv/dev/bin/framedepth | Statements Executed | 20741429 | Total Time | 54.9505060001434 seconds |
| Calls | P | F | Exclusive Time |
Inclusive Time |
Subroutine | |
|---|---|---|---|---|---|---|
| 984544 | 1 | 2 | 4.83287 | 4.83287 | main:: | gnomrev_c (xsub) |
| 984543 | 1 | 1 | 4.47873 | 9.31159 | main:: | gnomrev |
| 984543 | 1 | 1 | 4.15710 | 7.88087 | main:: | gnomfwd |
| 984544 | 1 | 2 | 3.72377 | 3.72377 | main:: | gnomfwd_c (xsub) |
| 0 | 0 | 0 | 0 | 0 | main:: | BEGIN |
| 0 | 0 | 0 | 0 | 0 | main:: | gnomfwd_pure |
| 0 | 0 | 0 | 0 | 0 | main:: | gnomrev_pure |
| 0 | 0 | 0 | 0 | 0 | main:: | picker |
| Line | Stmts. | Exclusive Time | Avg. | Code |
|---|---|---|---|---|
| 1 | 1 | 1.0e-6 | 1.0e-6 | #! /usr/bin/env perl |
| 2 | ||||
| 3 | 3 | 0.00119 | 0.00040 | use strict; # spent 12µs making 1 call to strict::import |
| 4 | 3 | 0.00222 | 0.00074 | use warnings; # spent 31µs making 1 call to warnings::import |
| 5 | ||||
| 6 | use WISE::Env ( | |||
| 7 | cfg => '/wise/base/deliv/dev/lib/perl', # spent 128ms making 1 call to WISE::Env::import | |||
| 8 | iam => 'FramePlot', | |||
| 9 | version=> '$Id: frameplot 5619 2009-08-21 20:30:50Z tim $', | |||
| 10 | params => [ | |||
| 11 | 'frame_table,frt,ft,t,f: file = $unnamed $required', | |||
| 12 | ' # Table of frames to plot', | |||
| 13 | 'coord_sys,coord,cs: list of 1 names (equ,ecl,gal)'. | |||
| 14 | ' = ecl', | |||
| 15 | ' # Coordinate system of plot', | |||
| 16 | 'zoom_to_fit,tofit,fit: switch', | |||
| 17 | ' # Zoom plot to just fit all frames', | |||
| 18 | 'window,win,w: list of 4 real', | |||
| 19 | ' # Sky window to plot: minlon,maxlon,minlat,maxlat', | |||
| 20 | ' # Absense means all-sky unless -fit specified', | |||
| 21 | 'pick,interact,int,i: switch', | |||
| 22 | ' # Enable interactive point picking (/xs device only)', | |||
| 23 | 'smooth,sm,s: switch', | |||
| 24 | ' # Smooth bin coverage estimes by averaging with neighbors', | |||
| 25 | 'res: list of 1 named (l,lo,low,m,med,h,hi,high) = hi', | |||
| 26 | ' # Desired resolution', | |||
| 27 | 'where: str', | |||
| 28 | ' # ADVANCED: Row selection criteria as pseudo-Perl', | |||
| 29 | 'getfix: switch', | |||
| 30 | ' # Run the standard frame index query and write the given '. | |||
| 31 | '-frame_file: SLOW!!!', | |||
| 32 | ' # Override/augment frmquery params with "... -:getfix \'...\' ..."', | |||
| 33 | 'plot_device,pdev,pd: str = /xs', | |||
| 34 | ' # PGPLOT device specification', | |||
| 35 | 'verbose,v: switch', | |||
| 36 | 'debug,dbg: str = 0', | |||
| 37 | 'plot_options,plotopts,popts: map (@defaults)'. | |||
| 38 | ' = zerocen=1,pole=0,ttl=0,corder=-1,'. | |||
| 39 | 'pad=0,binfill=undef,papx=0,binmark=-4,asp=0,'. | |||
| 40 | 'cmap=4,8,10,12,14,16,20,50,100,1000,', | |||
| 41 | ' # Plot customizations and options', | |||
| 42 | ' # cmap is the low end of each bin except the '. | |||
| 43 | 'lowest (depth>=1) bin', | |||
| 44 | 'options,opts: map (@defaults)'. | |||
| 45 | ' = subdiv=undef,binsz=undef,tblwcs=1,'. | |||
| 46 | 'scancol=scan,framecol=frame,inland=60,'. | |||
| 47 | 'wcscols=elon,elat,pixscl=2.75,npix=1016,'. | |||
| 48 | 'ixnsides=10', | |||
| 49 | ' # General options', | |||
| 50 | ], | |||
| 51 | import => [qw/%pvals $pars $err $warn/], | |||
| 52 | 3 | 0.00245 | 0.00082 | ); |
| 53 | ||||
| 54 | 3 | 5.9e-5 | 2.0e-5 | use File::Basename; # spent 86µs making 1 call to Exporter::import |
| 55 | 3 | 0.00053 | 0.00018 | use POSIX qw/asin/; # spent 2.29ms making 1 call to POSIX::import |
| 56 | ||||
| 57 | 3 | 0.00104 | 0.00035 | use WISE::CoUtils ('cconv','$R2D'); # spent 188µs making 1 call to Exporter::import |
| 58 | 3 | 0.00097 | 0.00032 | use WISE::Wrap; # spent 27µs making 1 call to Exporter::Lite::import |
| 59 | 3 | 0.00101 | 0.00034 | use WISE::PgAitoff 'make_aitoff'; # spent 66µs making 1 call to Exporter::import |
| 60 | 3 | 0.00137 | 0.00046 | use WISE::FITSIO::Utils; # spent 5µs making 1 call to import |
| 61 | 3 | 0.00074 | 0.00025 | use WISE::CHealPix; # spent 38µs making 1 call to Exporter::import |
| 62 | ||||
| 63 | #my @inline_dir; | |||
| 64 | #BEGIN { @inline_dir = (DIRECTORY => | |||
| 65 | # ("/wise/base/deliv/dev/lib/perl" !~ /<:/ ? "/wise/base/deliv/dev/lib/perl" : ".")."/_Inline"); } | |||
| 66 | ||||
| 67 | use Inline (C => 'DATA', # spent 134µs making 1 call to Inline::import | |||
| 68 | #AUTO_INCLUDE => '#include "rice.h"', | |||
| 69 | #INC => "-I/wise/base/deliv/dev/lib/perl/../../include", | |||
| 70 | #LIBS => "-L/wise/base/deliv/dev/lib/perl/..", | |||
| 71 | CLEAN_AFTER_BUILD => 0, | |||
| 72 | CLEAN_BUILD_AREA => 0, | |||
| 73 | #@inline_dir, | |||
| 74 | 3 | 0.00534 | 0.00178 | ); |
| 75 | ||||
| 76 | 1 | 2.0e-6 | 2.0e-6 | my @colrs = (2,5,6,8); # red,cyan,red-purple,orange |
| 77 | ||||
| 78 | # Appropriate for the number of colors expected; relies on internal | |||
| 79 | # knowledge of the PgAitoff module. | |||
| 80 | 1 | 4.0e-6 | 4.0e-6 | my @colptiles = (0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1); |
| 81 | 1 | 3.0e-6 | 3.0e-6 | my $cs = $pvals{coord_sys}[0]; |
| 82 | ||||
| 83 | 1 | 2.0e-6 | 2.0e-6 | $pvals{verbose} ||= $pvals{debug}; |
| 84 | ||||
| 85 | 1 | 5.0e-6 | 5.0e-6 | if($pvals{res}[0] =~ /^l/) { |
| 86 | 1 | 3.0e-6 | 3.0e-6 | $pvals{options}{binsz} //= 5600; |
| 87 | 1 | 1.0e-6 | 1.0e-6 | $pvals{options}{subdiv} //= 1; |
| 88 | 1 | 3.0e-6 | 3.0e-6 | $pvals{plot_options}{binfill} //= $cs eq 'ecl' ? 0.67 : 1; |
| 89 | } elsif($pvals{res}[0] =~ /^m/) { | |||
| 90 | $pvals{options}{binsz} //= 2800; | |||
| 91 | $pvals{options}{subdiv} //= 2; | |||
| 92 | $pvals{plot_options}{binfill} //= $cs eq 'ecl' ? 0.67 : 1; | |||
| 93 | } else { # high | |||
| 94 | $pvals{options}{binsz} //= 1400; | |||
| 95 | $pvals{options}{subdiv} //= 3; | |||
| 96 | $pvals{plot_options}{binfill} //= $cs eq 'ecl' ? 0.67 : 1; | |||
| 97 | } | |||
| 98 | ||||
| 99 | 1 | 1.4e-5 | 1.4e-5 | my $cdelt = $pvals{options}{pixscl}/3600; |
| 100 | 1 | 1.0e-6 | 1.0e-6 | my $npix = $pvals{options}{npix}; |
| 101 | 1 | 2.0e-6 | 2.0e-6 | my $cenpix= $npix/2 + 0.5; |
| 102 | 1 | 1.0e-6 | 1.0e-6 | my $width = $npix*$cdelt*3600; |
| 103 | 1 | 2.0e-6 | 2.0e-6 | my $inland= $pvals{options}{inland}; |
| 104 | 1 | 1.0e-6 | 1.0e-6 | my $inlandpix = $inland/($cdelt*3600); |
| 105 | 1 | 1.0e-6 | 1.0e-6 | my $pole = $pvals{options}{zerocen}; # not used right now |
| 106 | 1 | 0 | 0 | my $ncmap = 10; |
| 107 | 1 | 2.0e-6 | 2.0e-6 | my $cmap = $pvals{plot_options}{cmap} || []; |
| 108 | 1 | 1.0e-6 | 1.0e-6 | die "$err: Plot option cmap does't have 10 elements.\n" |
| 109 | if @$cmap != $ncmap; | |||
| 110 | 1 | 1.0e-6 | 1.0e-6 | die "$err: Can't do interactive except with the /xs pgplot device.\n" |
| 111 | if $pvals{pick} && $pvals{plot_device}!~m|/xs(erver)?$|; | |||
| 112 | 1 | 2.0e-6 | 2.0e-6 | my ($lonmin,$lonmax,$latmin,$latmax) = @{$pvals{window} || []}; |
| 113 | 1 | 2.0e-6 | 2.0e-6 | my $fit = $pvals{zoom_to_fit}; |
| 114 | 1 | 1.0e-6 | 1.0e-6 | my $zoom = defined $lonmin || $fit; |
| 115 | 1 | 1.0e-6 | 1.0e-6 | my $ixnsides = $pvals{options}{ixnsides} // 10; # Nsides to use for index binning |
| 116 | ||||
| 117 | 1 | 2.0e-6 | 2.0e-6 | my $tblfile = $pvals{frame_table} || '-'; |
| 118 | ||||
| 119 | 1 | 1.0e-6 | 1.0e-6 | if($pvals{getfix}) { |
| 120 | print "Running frame index query (be very patient) ... ".gmtime()." GMT\n"; | |||
| 121 | WISE::Wrap::getfix({ | |||
| 122 | f => {pstat=>0,score=>5}, | |||
| 123 | opts => {rowlimit=>0}, | |||
| 124 | }, | |||
| 125 | { | |||
| 126 | defs=>$pars, | |||
| 127 | noread=>1, | |||
| 128 | file=>$tblfile, | |||
| 129 | } | |||
| 130 | ) | |||
| 131 | or die "$err: Unable to search FIX.\n"; | |||
| 132 | } | |||
| 133 | ||||
| 134 | 1 | 7.0e-6 | 7.0e-6 | my @wcscols = map {lc} (@{$pvals{options}{wcscols}}); |
| 135 | 1 | 4.0e-6 | 4.0e-6 | my %wcscols = map { ($_=>1) } @wcscols; |
| 136 | 1 | 1.0e-6 | 1.0e-6 | my @eclcols = (qw/elon elat/); |
| 137 | 1 | 2.0e-6 | 2.0e-6 | my $ecl = $wcscols[0] =~ /^e/i; |
| 138 | 1 | 6.0e-6 | 6.0e-6 | my @frcols = map {lc} ($pvals{options}{scancol},$pvals{options}{framecol}) |
| 139 | if $pvals{options}{framecol} && $pvals{options}{scancol}; | |||
| 140 | #my @cols = (@frcols,@wcscols); | |||
| 141 | #print "Reading columns @cols from '$tblfile' ... ".gmtime()." GMT\n"; | |||
| 142 | 1 | 3.0e-5 | 3.0e-5 | print "Reading from table '$tblfile' ... ".gmtime()." GMT\n"; |
| 143 | 1 | 1.7e-5 | 1.7e-5 | my $tbl = WISE::IPACTbl->new($tblfile,{fast=>1}) # cols=>\@cols, # spent 8.76s making 1 call to WISE::IPACTbl::new |
| 144 | or die; | |||
| 145 | 1 | 3.7e-5 | 3.7e-5 | die "$err: Required columns @wcscols abesnt from '$tblfile'.\n" # spent 33µs making 1 call to WISE::IPACTbl::columns_read |
| 146 | if @wcscols && grep($wcscols{$_},$tbl->cols()) != @wcscols; | |||
| 147 | 1 | 9.0e-6 | 9.0e-6 | my $frs = $tbl->data(); # spent 102µs making 1 call to WISE::IPACTbl::data |
| 148 | 1 | 8.0e-6 | 8.0e-6 | my $nr = $tbl->nrows(); # spent 9µs making 1 call to WISE::IPACTbl::rows |
| 149 | 1 | 1.0e-6 | 1.0e-6 | die "$err: No frame info read from '$tblfile'.\n" |
| 150 | if $nr <= 0; | |||
| 151 | 1 | 4.5e-5 | 4.5e-5 | print "... read $nr framests.\n"; |
| 152 | ||||
| 153 | 1 | 1.2e-5 | 1.2e-5 | $tbl->addcol({name=>'_mch',type=>'i',blank=>'null'}); # spent 246µs making 1 call to WISE::IPACTbl::addcol |
| 154 | ||||
| 155 | 1 | 2.8e-5 | 2.8e-5 | my $nsides = WISE::CHealPix::side2nsides_deg($pvals{options}{binsz}/3600); # spent 23µs making 1 call to WISE::CHealPix::side2nsides_deg |
| 156 | 1 | 2.6e-5 | 2.6e-5 | my $binsz = WISE::CHealPix::c2cdist_deg($nsides)*3600; # spent 22µs making 1 call to WISE::CHealPix::c2cdist_deg |
| 157 | ||||
| 158 | 1 | 0 | 0 | my %bins; |
| 159 | 1 | 0 | 0 | my %index; |
| 160 | ||||
| 161 | 1 | 1.0e-6 | 1.0e-6 | my $n = 0; |
| 162 | 1 | 2.0e-6 | 2.0e-6 | my $nsteps = $pvals{options}{subdiv}; |
| 163 | 1 | 1.0e-6 | 1.0e-6 | my $d = $npix/$nsteps; |
| 164 | 1 | 0 | 0 | my $nfr = 0; |
| 165 | 1 | 1.0e-6 | 1.0e-6 | my $nbins = 0; |
| 166 | 1 | 1.0e-6 | 1.0e-6 | my ($frx0,$fry0) = (0.5,0.5); |
| 167 | ||||
| 168 | 1 | 0 | 0 | my $wheresub; |
| 169 | 1 | 1.0e-6 | 1.0e-6 | if($pvals{where}) { |
| 170 | # Replace $[colname] with a column value | |||
| 171 | (my $where = $pvals{where}) =~ s/\$\[\s*([^]]+?)\s*\]/\$rows->{$1}[\$row]/g; | |||
| 172 | my $code = "sub {my \$rows=shift; my \$row=shift; $where;}"; | |||
| 173 | $wheresub = WISE::Utils::safe_eval($code, | |||
| 174 | {share_from=> | |||
| 175 | ['main' => [qw/&rad &hypot/], | |||
| 176 | 'WISE::CoUtils' => [qw/&angdist/], | |||
| 177 | 'WISE::CoUtilsXS'=> [], | |||
| 178 | ], | |||
| 179 | permit=>[qw/grepstart grepwhile/]}); | |||
| 180 | die "$err: Couldn't eval where code '$code' from '$where'.\n$@" | |||
| 181 | if $@; | |||
| 182 | print "Applying row selection routine '$code' ...\n"; | |||
| 183 | } | |||
| 184 | ||||
| 185 | 1 | 6.3e-5 | 6.3e-5 | print "Binning $nr framesets ". |
| 186 | "(step=$d pixels, nsteps=$nsteps, binsz=$binsz\", nsides=$nsides) ... ".gmtime()." GMT\n"; | |||
| 187 | ||||
| 188 | 1 | 3.0e-6 | 3.0e-6 | for my $i (0..$nr-1) { |
| 189 | 984543 | 4.30170 | 4.4e-6 | my ($elon0,$elat0) = map {$frs->{$_}[$i]} @wcscols; |
| 190 | 984543 | 0.54291 | 5.5e-7 | next if ! defined $elon0; |
| 191 | 984543 | 0.87091 | 8.8e-7 | $frs->{_mch}[$i] = 0; |
| 192 | 984543 | 4.41294 | 4.5e-6 | my $ixbin = WISE::CHealPix::ang2pix_ring($ixnsides,$elon0,$elat0); # spent 12.1s making 984543 calls to WISE::CHealPix::ang2pix_ring, avg 12µs/call |
| 193 | 984543 | 1.21768 | 1.2e-6 | push @{$index{$ixbin}},$i; |
| 194 | 984543 | 0.52047 | 5.3e-7 | next if $wheresub && ! $wheresub->($frs,$i); |
| 195 | 984543 | 0.74166 | 7.5e-7 | $frs->{_mch}[$i] = 1; |
| 196 | 984543 | 0.42029 | 4.3e-7 | my %frbins; |
| 197 | 984543 | 1.67469 | 1.7e-6 | for my $ix (0..$nsteps-1) { |
| 198 | 984543 | 1.16388 | 1.2e-6 | my $x = $frx0+$ix*$d + $d/2; |
| 199 | 984543 | 2.12439 | 2.2e-6 | for my $iy (0..$nsteps-1) { |
| 200 | 984543 | 0.80727 | 8.2e-7 | my $y = $fry0+$iy*$d + $d/2; |
| 201 | 984543 | 4.66981 | 4.7e-6 | my ($elon,$elat) = gnomrev($x,$y,$elon0,$elat0); # spent 9.31s making 984543 calls to main::gnomrev, avg 9µs/call |
| 202 | 984543 | 3.96330 | 4.0e-6 | my $bin = WISE::CHealPix::ang2pix_ring($nsides,$elon,$elat); # spent 10.3s making 984543 calls to WISE::CHealPix::ang2pix_ring, avg 10µs/call |
| 203 | 984543 | 3.02310 | 3.1e-6 | if(! $frbins{$bin}++) { |
| 204 | 984543 | 4.32619 | 4.4e-6 | ($elon,$elat) = WISE::CHealPix::pix2ang_ring($nsides,$bin); # spent 11.6s making 984543 calls to WISE::CHealPix::pix2ang_ring, avg 12µs/call |
| 205 | 984543 | 4.34683 | 4.4e-6 | my ($xp,$yp) = gnomfwd($elon,$elat,$elon0,$elat0); # spent 7.88s making 984543 calls to main::gnomfwd, avg 8µs/call |
| 206 | 984543 | 1.50665 | 1.5e-6 | if($xp>$inlandpix+0.5 && $xp<1016.5-$inlandpix && |
| 207 | $yp>$inlandpix+0.5 && $yp<1016.5-$inlandpix) { | |||
| 208 | ++$bins{$bin}; | |||
| 209 | } | |||
| 210 | } | |||
| 211 | } | |||
| 212 | } | |||
| 213 | 984543 | 1.38877 | 1.4e-6 | ++$nfr; |
| 214 | } | |||
| 215 | ||||
| 216 | 1 | 4.0e-6 | 4.0e-6 | $nbins = keys %bins; |
| 217 | ||||
| 218 | 1 | 3.0e-6 | 3.0e-6 | if($pvals{smooth}) { |
| 219 | print "Smoothing $nbins coverage bins ... ".gmtime()." GMT\n"; | |||
| 220 | my %newbins; | |||
| 221 | for my $bin (keys %bins) { | |||
| 222 | my @bins = WISE::CHealPix::neighbors_ring($nsides,$bin); | |||
| 223 | my $ave = 0; | |||
| 224 | my $n = 0; | |||
| 225 | for ($bin,@bins) { | |||
| 226 | $ave += $bins{$_}//0; | |||
| 227 | ++$n; | |||
| 228 | } | |||
| 229 | $newbins{$bin} = int($ave/$n + 0.5); | |||
| 230 | } | |||
| 231 | %bins = %newbins; | |||
| 232 | } | |||
| 233 | ||||
| 234 | 1 | 3.2e-5 | 3.2e-5 | print "Getting bin centers for $nfr frames ... ".gmtime()." GMT\n"; |
| 235 | ||||
| 236 | 1 | 1.0e-6 | 1.0e-6 | my (@x,@y,@z); |
| 237 | 1 | 0.00510 | 0.00510 | for my $bin (keys %bins) { |
| 238 | 14018 | 0.06227 | 4.4e-6 | my ($elon,$elat) = WISE::CHealPix::pix2ang_ring($nsides,$bin); # spent 174ms making 14018 calls to WISE::CHealPix::pix2ang_ring, avg 12µs/call |
| 239 | 14018 | 0.00918 | 6.6e-7 | push @x,$elon; |
| 240 | 14018 | 0.00766 | 5.5e-7 | push @y,$elat; |
| 241 | 14018 | 0.02979 | 2.1e-6 | push @z,$bins{$bin}; |
| 242 | ||||
| 243 | } | |||
| 244 | ||||
| 245 | 1 | 2.1e-5 | 2.1e-5 | print "Coverage histogram:\n"; |
| 246 | ||||
| 247 | 1 | 1.0e-6 | 1.0e-6 | my @color_key; |
| 248 | 1 | 1.0e-6 | 1.0e-6 | my $cov = 1; |
| 249 | 1 | 1.7e-5 | 1.7e-5 | for my $i (0..$#$cmap,$#$cmap+1) { |
| 250 | 11 | 0.02668 | 0.00243 | my $n = grep { $_ >= $cov } values %bins; |
| 251 | 11 | 3.3e-5 | 3.0e-6 | my $area = $n*($binsz/3600)**2; |
| 252 | 11 | 0.00026 | 2.3e-5 | printf " %8.2f sq. deg.s >= $cov coverages\n",$area,$cov; |
| 253 | 11 | 5.9e-5 | 5.4e-6 | $color_key[$i+1] = sprintf("%4d".($zoom?"":" (%4.0f)"),$cov,$area); |
| 254 | 11 | 2.8e-5 | 2.5e-6 | $cov = $cmap->[$i]; |
| 255 | } | |||
| 256 | ||||
| 257 | 1 | 1.0e-6 | 1.0e-6 | my ($min,$max); |
| 258 | 1 | 0.00400 | 0.00400 | ($min,$max) = (sort {$a<=>$b} @z)[0,$#z] if @z; |
| 259 | ||||
| 260 | 1 | 5.6e-5 | 5.6e-5 | print "Colorizing from $min to $max ... ".gmtime()." GMT\n"; |
| 261 | ||||
| 262 | # Map to 11 colors by specified range | |||
| 263 | 1 | 1.0e-6 | 1.0e-6 | my %color_map; |
| 264 | 1 | 6.0e-6 | 6.0e-6 | for my $depth ($min..$max) { |
| 265 | 589 | 0.00024 | 4.1e-7 | my $c; |
| 266 | 589 | 0.00070 | 1.2e-6 | for my $i (0..$#$cmap) { |
| 267 | 5665 | 0.00504 | 8.9e-7 | if($depth < $cmap->[$i]) { |
| 268 | 589 | 0.00027 | 4.5e-7 | $c = $i+1; |
| 269 | 589 | 0.00040 | 6.8e-7 | last; |
| 270 | } | |||
| 271 | } | |||
| 272 | 589 | 0.00026 | 4.3e-7 | $c ||= $ncmap + 1; |
| 273 | 589 | 0.00110 | 1.9e-6 | $color_map{$depth} = $c; |
| 274 | 589 | 0.00067 | 1.1e-6 | print " $depth -> $c\n" if $pvals{debug}; |
| 275 | } | |||
| 276 | ||||
| 277 | 1 | 1.0e-6 | 1.0e-6 | my ($x0,$x1,$x0z,$x1z,$y0,$y1); |
| 278 | 1 | 1.2e-5 | 1.2e-5 | if($cs ne 'ecl' || $zoom) { |
| 279 | print "Converting to coordinate '$cs' ... ".gmtime()." GMT\n" | |||
| 280 | if $cs ne 'ecl';; | |||
| 281 | for (0..$#x) { | |||
| 282 | ($x[$_],$y[$_]) = WISE::CoUtils::cconv('ecl',$x[$_],$y[$_],$cs) | |||
| 283 | if $cs ne 'ecl'; | |||
| 284 | if($zoom) { | |||
| 285 | my $xz = $x[$_]>180 ? $x[$_] - 360 : $x[$_]; | |||
| 286 | $x0 = ! defined $x0 || $x[$_]<$x0 ? $x[$_] : $x0; | |||
| 287 | $x1 = ! defined $x1 || $x[$_]>$x1 ? $x[$_] : $x1; | |||
| 288 | $x0z = ! defined $x0z || $xz <$x0z ? $xz : $x0z; | |||
| 289 | $x1z = ! defined $x1z || $xz >$x1z ? $xz : $x1z; | |||
| 290 | $y0 = ! defined $y0 || $y[$_]<$y0 ? $y[$_] : $y0; | |||
| 291 | $y1 = ! defined $y1 || $y[$_]>$y1 ? $y[$_] : $y1; | |||
| 292 | } | |||
| 293 | } | |||
| 294 | } | |||
| 295 | ||||
| 296 | 1 | 1.0e-6 | 1.0e-6 | die "$err: No data.\n" if ! @x; |
| 297 | ||||
| 298 | ||||
| 299 | # Cross 0/360? | |||
| 300 | 1 | 4.0e-6 | 4.0e-6 | my $cen = $pvals{plot_options}{zerocen} || ($zoom && $x1-$x0 > 180); |
| 301 | ||||
| 302 | 1 | 1.0e-6 | 1.0e-6 | if($fit) { |
| 303 | $lonmin = $cen ? $x0z : $x0; | |||
| 304 | $lonmax = $cen ? $x1z : $x1; | |||
| 305 | $latmin = $y0; | |||
| 306 | $latmax = $y1; | |||
| 307 | } | |||
| 308 | ||||
| 309 | 1 | 1.0e-6 | 1.0e-6 | if($fit && $lonmax - $lonmin >= 180) { |
| 310 | # Stick with all-sky | |||
| 311 | $lonmin = undef; | |||
| 312 | $zoom = 0; | |||
| 313 | } | |||
| 314 | ||||
| 315 | 1 | 1.9e-5 | 1.9e-5 | print "\nPlotting to $pvals{plot_device} ... ".gmtime()." GMT\n"; |
| 316 | ||||
| 317 | 1 | 1.0e-6 | 1.0e-6 | print "\nInteractive keys: any=position, q=quit, f=overlapping frames, d=depth\n\n" |
| 318 | if $pvals{pick}; | |||
| 319 | ||||
| 320 | 1 | 9.0e-5 | 9.0e-5 | make_aitoff({ # spent 772ms making 1 call to WISE::PgAitoff::make_aitoff |
| 321 | x => \@x, # Bin/marker lon | |||
| 322 | 'y' => \@y, # Bin/marker lat | |||
| 323 | z => \@z, # Bin/marker colors | |||
| 324 | lon0 => $pvals{plot_options}{zerocen} ? 0 : 180, | |||
| 325 | x0 => $lonmin, # Coord. grid zoom (if defined) | |||
| 326 | x1 => $lonmax, | |||
| 327 | y0 => $latmin, | |||
| 328 | y1 => $latmax, | |||
| 329 | binsz => $binsz/3600, # Bin dimension | |||
| 330 | binfill => $pvals{plot_options}{binfill}, # Bin fill factor | |||
| 331 | spec => 1, # Pseudo thermal/spectral color map | |||
| 332 | nolow => 1, # No "too low" bin in color map | |||
| 333 | keyttl => 'Depth'.($zoom?'':' (sq. deg.s)'), # Z-axis title | |||
| 334 | keylbls => \@color_key, # Color key; only used for depth display | |||
| 335 | colmap => \%color_map, | |||
| 336 | corder => $pvals{plot_options}{corder}, # Plot colors in reverse depth-coverage order | |||
| 337 | coord => $cs, | |||
| 338 | pole => $pole, | |||
| 339 | title => $pvals{plot_options}{ttl} || "Exposure Depth (frames)", | |||
| 340 | subtitle => "NFrames = $nfr, File = $tblfile", | |||
| 341 | pick => $pvals{pick} && \&picker, | |||
| 342 | papxwidth=> $pvals{plot_options}{papx}, | |||
| 343 | binmark => $pvals{plot_options}{binmark}, | |||
| 344 | plotasp => $pvals{plot_options}{asp}, | |||
| 345 | pdev => $pvals{plot_device}, | |||
| 346 | verbose => $pvals{debug}, | |||
| 347 | debug => $pvals{debug}, | |||
| 348 | }); | |||
| 349 | ||||
| 350 | 1 | 4.2e-5 | 4.2e-5 | print "\nDone. ".gmtime()." GMT\n"; |
| 351 | ||||
| 352 | 1 | 7.0e-5 | 7.0e-5 | exit 0; |
| 353 | ||||
| 354 | ||||
| 355 | # ========== | |||
| 356 | # Subs | |||
| 357 | # ========== | |||
| 358 | ||||
| 359 | # Look up input frames (or bins) by target coord lon,lat. | |||
| 360 | # Callback for make_aitoff. | |||
| 361 | sub picker { | |||
| 362 | my ($c,$x,$y) = @_; | |||
| 363 | return 0 if $c =~ /^[q]$/i; # Done | |||
| 364 | next if ! defined $x; | |||
| 365 | my $rad = ($width*sqrt(2)+$inland)/3600; | |||
| 366 | my ($ra,$dec) = WISE::CoUtils::cconv($cs,$x,$y,'equ'); | |||
| 367 | my ($elon,$elat) = WISE::CoUtils::cconv('equ',$ra,$dec,'ecl'); | |||
| 368 | my ($glon,$glat) = WISE::CoUtils::cconv('equ',$ra,$dec,'gal'); | |||
| 369 | printf "\n(RA, Dec) = %10.6f,%10.6f\n",$ra,$dec; | |||
| 370 | printf "(Glon,Glat) = %10.6f,%10.6f\n",$glon,$glat; | |||
| 371 | printf "(Elon,Elat) = %10.6f,%10.6f\n",$elon,$elat; | |||
| 372 | my $bin = WISE::CHealPix::ang2pix_ring($nsides,$elon,$elat); | |||
| 373 | print "Approximate, binned depth (frames) = ".($bins{$bin}//0)."\n"; | |||
| 374 | # Look-up overlapping frames | |||
| 375 | if($c =~ /[fd]/i) { | |||
| 376 | my $mch = 0; | |||
| 377 | my @mchs; | |||
| 378 | print "Searching for overlapping frames ...\n"; | |||
| 379 | my $ixbin = WISE::CHealPix::ang2pix_ring($ixnsides,$elon,$elat); | |||
| 380 | my @ixbins = WISE::CHealPix::pix_withinr_ring($ixnsides,$ixbin,4200/3600); | |||
| 381 | my @i; | |||
| 382 | push @i, @{$index{$_} || []} for @ixbins; | |||
| 383 | for my $i (sort {$a<=>$b} @i) { | |||
| 384 | my ($elon0,$elat0) = map {$frs->{$_}[$i]} @eclcols; | |||
| 385 | next if ! defined $elon0; | |||
| 386 | next if abs($elat-$elat0) > $rad; | |||
| 387 | next if WISE::CoUtils::angdist($elon,$elat,$elon0,$elat0) > $rad; | |||
| 388 | my ($x,$y) = gnomfwd($elon, $elat, $elon0, $elat0); | |||
| 389 | my $xin= $x<$cenpix ? $x - 0.5 : $npix - $x + 0.5; | |||
| 390 | my $yin= $y<$cenpix ? $y - 0.5 : $npix - $y + 0.5; | |||
| 391 | my $edgedist= $xin<$yin ? -$xin*$cdelt*3600 : -$yin*$cdelt*3600; | |||
| 392 | if($edgedist <= -($pvals{options}{inland}//0)) { | |||
| 393 | my ($scan,$frame) = map {$frs->{$_}[$i]} @frcols; | |||
| 394 | ++$mch if $frs->{_mch}[$i]; | |||
| 395 | if($c eq 'f') { | |||
| 396 | # Short summary of frame info | |||
| 397 | printf " %4d: scan/frame=%s%03d, inland dist=%6.1f\", ". | |||
| 398 | "(elon,elat=%10.6f,%10.6f,x,y=%4.0f,%4.0f)\n", | |||
| 399 | $mch,$scan,$frame,-$edgedist,$elon0,$elat0,$x,$y; | |||
| 400 | } elsif($c eq 'F') { | |||
| 401 | # Save matched frame FIX info | |||
| 402 | push @mchs, { map { (lc($_)=>$frs->{lc $_}[$i]) } $tbl->names() }; | |||
| 403 | } | |||
| 404 | } | |||
| 405 | } | |||
| 406 | if(@mchs) { | |||
| 407 | # Display matched frame info | |||
| 408 | WISE::IPACTbl->new("-","w",{meta=>$tbl,data=>\@mchs})->data(\@mchs); | |||
| 409 | } | |||
| 410 | print "True depth (frames) = $mch\n"; | |||
| 411 | } | |||
| 412 | print "\nInteractive keys: <any>=position, q=quit, f/F=overlapping frames, d=depth\n"; | |||
| 413 | return 1; | |||
| 414 | } | |||
| 415 | ||||
| 416 | 984543 | 5.79067 | 5.9e-6 | # spent 7.88s (4.16+3.72) within main::gnomfwd which was called 984543 times, avg 8µs/call:
# 984543 times (4.16s+3.72s) at line 205, avg 8µs/call # spent 3.72s making 984543 calls to main::gnomfwd_c, avg 4µs/call |
| 417 | 984543 | 6.96502 | 7.1e-6 | # spent 9.31s (4.48+4.83) within main::gnomrev which was called 984543 times, avg 9µs/call:
# 984543 times (4.48s+4.83s) at line 201, avg 9µs/call # spent 4.83s making 984543 calls to main::gnomrev_c, avg 5µs/call |
| 418 | ||||
| 419 | #*gnomfwd = \&gnomfwd_pure; | |||
| 420 | #*gnomrev = \&gnomrev_pure; | |||
| 421 | ||||
| 422 | # Pure perl gnomonic projections | |||
| 423 | # COPIED and modified from CoUtils.pm to reduce arg count | |||
| 424 | sub gnomrev_pure { | |||
| 425 | my ($xi,$eta,$ra0,$dec0) = @_; | |||
| 426 | my $xp = ($xi - $cenpix)*$cdelt/$R2D; | |||
| 427 | my $yp = ($eta - $cenpix)*$cdelt/$R2D; | |||
| 428 | # delta = atan(sqrt($xp*$xp + $yp*$yp)) | |||
| 429 | my $rsq = $xp*$xp + $yp*$yp; | |||
| 430 | my $r = sqrt($rsq); | |||
| 431 | my $cosdelta = $rsq > 1e-32 ? 1. / sqrt($rsq + 1) : 1.0; | |||
| 432 | my $sindelta = sqrt(1. - $cosdelta*$cosdelta); | |||
| 433 | # beta = atan2(-$x,$y) | |||
| 434 | my $cosbeta = $r>0 ? -$yp/$r : 1; | |||
| 435 | my $sinbeta = $r>0 ? -$xp/$r : 0; | |||
| 436 | my $cosdec0 = cos($dec0/$R2D); | |||
| 437 | my $sindec0 = sin($dec0/$R2D); | |||
| 438 | my $xx = $sindec0*$sindelta*$cosbeta + $cosdec0*$cosdelta; | |||
| 439 | my $yy = $sindelta*$sinbeta; | |||
| 440 | my $dec = asin($sindec0*$cosdelta - $cosdec0*$sindelta*$cosbeta)*$R2D; | |||
| 441 | my $ra = $ra0 + atan2($yy,$xx)*$R2D; | |||
| 442 | if($ra >= 360) { $ra -= 360 ; } | |||
| 443 | if($ra < 0) { $ra += 360 ; } | |||
| 444 | #printf("--- rev: returning %f,%f.\n",$ra,$dec); | |||
| 445 | return ($ra,$dec); | |||
| 446 | } | |||
| 447 | ||||
| 448 | sub gnomfwd_pure { | |||
| 449 | my ($ra, $dec, $ra0, $dec0) = @_; | |||
| 450 | my $cosdec0 = cos($dec0/$R2D); | |||
| 451 | my $sindec0 = sin($dec0/$R2D); | |||
| 452 | my $drar = ($ra0-$ra)/$R2D; | |||
| 453 | my $decr = $dec/$R2D; | |||
| 454 | my $a = cos($decr) * cos($drar); | |||
| 455 | my $f = 1.0/($sindec0 * sin($decr) + $a*$cosdec0); | |||
| 456 | my $x = $f * cos($decr) * sin($drar); | |||
| 457 | my $y = $f * ($cosdec0 * sin($decr) - $a*$sindec0); | |||
| 458 | my $xi = $x*$R2D/$cdelt + $cenpix; | |||
| 459 | my $eta = $y*$R2D/$cdelt + $cenpix; | |||
| 460 | return ($xi,$eta); | |||
| 461 | } | |||
| 462 | ||||
| 463 | ||||
| 464 | ########################## C subs | |||
| 465 | ||||
| 466 | __END__ | |||
| 467 | __C__ | |||
| 468 | ||||
| 469 | #define PI 3.1415926535897931 | |||
| 470 | #define R2D (180.0 / PI) | |||
| 471 | ||||
| 472 | void | |||
| 473 | gnomfwd_c(double ra, double dec, double ra0, double dec0, double cdelt, double cenpix) | |||
| 474 | { | |||
| 475 | Inline_Stack_Vars; | |||
| 476 | double cosdec0 = cos(dec0/R2D); | |||
| 477 | double sindec0 = sin(dec0/R2D); | |||
| 478 | double drar = (ra0-ra)/R2D; | |||
| 479 | double decr = dec/R2D; | |||
| 480 | double a = cos(decr) * cos(drar); | |||
| 481 | double f = 1.0/(sindec0 * sin(decr) + a*cosdec0); | |||
| 482 | double x = f * cos(decr) * sin(drar); | |||
| 483 | double y = f * (cosdec0 * sin(decr) - a*sindec0); | |||
| 484 | double xi = x*R2D/cdelt + cenpix; | |||
| 485 | double eta = y*R2D/cdelt + cenpix; | |||
| 486 | Inline_Stack_Reset; | |||
| 487 | Inline_Stack_Push(sv_2mortal(newSVnv(xi))); | |||
| 488 | Inline_Stack_Push(sv_2mortal(newSVnv(eta))); | |||
| 489 | Inline_Stack_Done; | |||
| 490 | } | |||
| 491 | ||||
| 492 | void | |||
| 493 | gnomrev_c(double xi, double eta, double ra0, double dec0, double cdelt, double cenpix) | |||
| 494 | { | |||
| 495 | Inline_Stack_Vars; | |||
| 496 | double xp = (xi - cenpix)*cdelt/R2D; | |||
| 497 | double yp = (eta - cenpix)*cdelt/R2D; | |||
| 498 | double rsq = xp*xp + yp*yp; | |||
| 499 | double r = sqrt(rsq); | |||
| 500 | double cosdelta = rsq > 1e-32 ? 1. / sqrt(rsq + 1) : 1.0; | |||
| 501 | double sindelta = sqrt(1. - cosdelta*cosdelta); | |||
| 502 | double cosbeta = r>0 ? -yp/r : 1; | |||
| 503 | double sinbeta = r>0 ? -xp/r : 0; | |||
| 504 | double cosdec0 = cos(dec0/R2D); | |||
| 505 | double sindec0 = sin(dec0/R2D); | |||
| 506 | double xx = sindec0*sindelta*cosbeta + cosdec0*cosdelta; | |||
| 507 | double yy = sindelta*sinbeta; | |||
| 508 | double dec = asin(sindec0*cosdelta - cosdec0*sindelta*cosbeta)*R2D; | |||
| 509 | double ra = ra0 + atan2(yy,xx)*R2D; | |||
| 510 | if(ra >= 360) { ra -= 360 ; } | |||
| 511 | if(ra < 0) { ra += 360 ; } | |||
| 512 | Inline_Stack_Reset; | |||
| 513 | Inline_Stack_Push(sv_2mortal(newSVnv(ra))); | |||
| 514 | Inline_Stack_Push(sv_2mortal(newSVnv(dec))); | |||
| 515 | Inline_Stack_Done; | |||
| 516 | } | |||
# spent 3.72s within main::gnomfwd_c which was called 984543 times, avg 4µs/call:
# 984543 times (3.72s+0) by main::gnomfwd at line 416 of /wise/base/deliv/dev/bin/framedepth, avg 4µs/call | ||||
# spent 4.83s within main::gnomrev_c which was called 984543 times, avg 5µs/call:
# 984543 times (4.83s+0) by main::gnomrev at line 417 of /wise/base/deliv/dev/bin/framedepth, avg 5µs/call |