PNe
#new enough version (PNe26) to make a fresh start,
but still rely
# heavily on the old PNe.sm code
(from PNe24, etc)
#read version file to figure out what pre-fix number
everything gets
data v.v
read { vnum 1 }
define v $(vnum[0])
#read iteration file to figure out the other pre-fix
number
data i.i
read { inum 1 }
define itn $(inum[0])
#advance the iteration number for the next instance
of this program
!rm i.i
write i.i $($itn + 1)
#NO FRILLS VERSION, STILL, NO CONTAM, PROPORTION,
ANYTHING
#SMOOTHING 9X9 - NO EXCEPTIONS (means r = 4)
#read in PNe datafile for feldmeier fields
data PNeF.dat
read { field 1.s num 2 RAh 3 RAm 4 RAs 5 DECd 6
DECam 7 DECas 8 m5007 9 notes 10.s }
#set RA/DEC both in degrees and range, too.
set DECf = DECd + (DECam/60) + (DECas/3600)
set RAf = 15 * (RAh + (RAm/60) + (RAs/3600))
#read in PNe datafile for aguerri fields
data PNeA.dat
read { field 1.s num 2 RAh 3 RAm 4 RAs 5 DECd 6
DECam 7 DECas 8 m5007 9 notes 10.s }
set DECa = DECd + (DECam/60) + (DECas/3600)
set RAa = 15 * (RAh + (RAm/60) + (RAs/3600))
#master setting of running parameters
set imOs = { 3.fits 3m.fits 4.fits 4m.fits 7.fits
7m.fits Core.fits \
Corem.fits Sub.fits Subm.fits
FCJ.fits FCJm.fits LPC.fits \
LPCm.fits }
set imO = { 0 1 2 3 4 5 6 7 8 9 10 11 12 13 }
set raO = { 1 3 5 7 9 11 13 15 25 35 45 55 65 }
set PNeNO = { 75 75 12 12 77 77 38 38 107 107 37 37
19 19 }
#temp for testing
# set imsO = { 3.fits 3m.fits 4.fits 4m.fits 7.fits
7m.fits Core.fits Corem.fits Sub.fits Subm.fits FCJ.fits FCJm.fits
LPC.fits LPCm.fits }
# set imO = { 0 }
# set imO =
{
1 2
3 4
5 6
7
8
9
10
11
12 13 }
#for reference 0
1 2
3 4
5
6 \
#
7
8
9
10
11
12 13 }
#PUT A NUMBER (0-13) UNDER IMAGES YOU WANT, NOTHING
UNDER IMAGES
# YOU DON'T WANT
# set raO = { 15 }
set raO = { 5 7 9 11 13 15 25 35 45 55 }
define smn 9
set PNeIM = { 3all$smn.fits 3mall$smn.fits
4all$smn.fits
\
4mall$smn.fits 7all$smn.fits
7mall$smn.fits Coreall$smn.fits \
Coremall$smn.fits Suball$smn.fits
Submall$smn.fits
\
FCJall$smn.fits FCJmall$smn.fits
LPCall$smn.fits
\
LPCmall$smn.fits }
#big images/ranges loops
foreach rang raO {
foreach imagN imO {
#set number of
PNe and master image to load
define imag
$(imOs[$imagN])
define PNeimag
$(PNeIM[$imagN])
define nPNe
$(PNeNO[$imagN])
#separating
feldmeier from aguerri more
if ('$imag' ==
'3.fits' || '$imag' == '3m.fits' || '$imag' == '4.fits' || '$imag' ==
'4m.fits' || '$imag' == '7.fits' || '$imag' == '7m.fits') {
set RA = RAf
set DEC = DECf
set rea = reaf
}
if ('$imag' ==
'Core.fits' || '$imag' == 'Corem.fits' || '$imag' == 'Sub.fits' ||
'$imag' == 'Subm.fits' || '$imag' == 'LPC.fits' || '$imag' ==
'LPCm.fits' || '$imag' == 'FCJ.fits' || '$imag' == 'FCJm.fits') {
set RA = RAa
set DEC = DECa
set rea = reaa
}
#read in
original image/fits stuff
define
file_type fits
image $imag
#prep fits/etc
data
#x-axis length
define NAXIS1
image
#y-axis length
define NAXIS2
image
#0 point in
RA/DEC
define CRVAL1
image
define CRVAL2
image
#0 point in
image x/y
define CRPIX1
image
define CRPIX2
image
#pixel scale
in dRA/dx and dDEC/dy
define CD1_1
image
define CD2_2
image
#scale range
from pixels into DEGREES
define rangeRD
(($rang)*ABS($CD1_1))
#now load in
the combined image, for mask-checking
image $PNeimag
#prep data
set bins = -5,
100, 0.5
set
DIMEN(wPNe) = 0
set
DIMEN(woPNe) = 0
set
DIMEN(rand) = 0
#files to save
box coordinates in and compare to later
#file for
boxes with 1+ PNe
define wfound
"$!v.$!itn.Im$!imag.R$!rang.wfound.dat"
#file for
boxes with 0 PNe
define zero
"$!v.$!itn.Im$!imag.R$!rang.0found.dat"
#file for
random boxes
define r
"$!v.$!itn.Im$!imag.R$!rang.Rfound.dat"
#prep vectors
to save coordinates
set dimen(saveRAw) = 0
set dimen(saveDECw) = 0
set
dimen(saveRAz) = 0
set
dimen(saveDECz) = 0
set
dimen(saveRAr) = 0
set
dimen(saveDECr) = 0
#keep track of
how many have been found, initialise to 0
define
foundPNe1 0
define
foundPNe0 0
define
foundPNeR 0
#run
with/without PNe boxes
echo calling
with and without PNe boxes
PNe10
#run random
boxes
echo calling
random boxes
PNeR
#run grapher
PNeG
#end of images loop
}
#end of ranges loop
}
pickgoodbox 4
#small macro to select a random box and return its
pixel values and RA,
# DEC values in the 2 passed
variable names
#choose random coordinates within $rang of the edge.
smoothing's done
# already so we don't care about
smoothing size
define xgood 0
define ygood 0
#make while loop that tries until a good coordinate
pair is found
while {$xgood == 0 || $ygood == 0} {
#generate randoms
define xtemp
$(int(($(RANDOM(1))*($NAXIS1-(2*$rang)))+$rang))
define ytemp
$(int(($(RANDOM(1))*($NAXIS2-(2*$rang)))+$rang))
#check if center is masked in mm
theoretical-generating
# smoothed image
if
($(image[$($xtemp+$NAXIS1),$ytemp]) > -700) {
#center is
unmasked, and we're ok to continue
#do first scan
over all pixels (from mm [+$NAXIS1]
#
theoretical-generating smoothed image) in box
#
to get number of masked pixels
define
nummasked 0
do x =
$($xtemp-$rang+$NAXIS1), $($xtemp+$rang+$NAXIS1) {
do y = $($ytemp-$rang), $($ytemp+$rang) {
if ($(image[$x,$y]) < -700) {
#this pixel's masked
define nummasked
$($nummasked+1)
}
}
}
#compute
masked fraction to use to adjust probability
define maskedf
$($nummasked/$($rang*2 + 1))
#we have a
winner, set xgood/ygood as them
define xgood
$xtemp
define ygood
$ytemp
#end of if the center is unmasked
}
#end of trying to find a good coordinate pair
}
#set inputs AND also set RA/DEC, just to be efficient
define $1 $xgood
define $2 $ygood
coo2wcs $1 $2 $3 $4
coo2wcs 4
#small macro to convert x,y coordinates into RA DEC
#coo2wcs x y RA DEC
define $4 $((($$2 - $CRPIX2 + 1)*$(ABS($CD2_2)) +
$CRVAL2))
define $3 $($(-($$1 - $CRPIX1 + 1
))*$(ABS($CD1_1)/cosd($$4)) + $CRVAL1)
wcs2coo 4
#small macro to convert RA DEC into xy
#wcs2coo RA DEC x y
PNe10
#boxes with and without PNe
#loop while either of these values is less than the
nPNe value
while {$foundPNe1 < $nPNe || $foundPNe0 <
$nPNe} {
#call pickbox with rx ry to get
random coordinates
pickgoodbox rx ry rRA rDEC
#this gives us a GOOD box, and
sets the maskedf for use later
#looking for PNe in a box
centered around the point that is 2*rang x 2*range
set delRA=$rRA-RA
set delDEC=$rDEC-DEC
set counter=( (abs(delRA)<=$rangeRD) &&
(abs(delDEC)<=$rangeRD) )? 1 : 0
if ($(sum(counter)) > 0) {
#found at
least one PN in this box
#check to be
sure we need it
if ($foundPNe1
< $nPNe) {
#we need it, let's use it
#save coords/pixels, iterate counter
set saveRAw = saveRAw concat $rRA
set saveDECw = saveDECw concat $rDEC
define foundPNe1 $($foundPNe1+1)
echo foundPNe1 = $foundPNe1
set dimen(statra)=$($(2*$rang+1)*$(2*$rang+1))
define idum 0
do x = -($rang), $rang {
do y = -($rang), $rang {
set statra[$idum]=$(image[$rx+$x,$ry+$y])
define idum $($idum+1)
}
}
#handles masked pixels
set statrafinal=statra if (statra>-900)
#concat this image onto the growing big file
set wPNe = wPNe concat statrafinal
#end of making
sure we need it
}
#end of found a PN
} ELSE {
#found 0 PN in
this box
#check that we
need it
if ($foundPNe0
< $nPNe) {
#we do! let's do it!
#save coords/pixels, iterate counter
set saveRAz = saveRAz concat $rRA
set saveDECz = saveDECz concat $rDEC
define foundPNe0 $($foundPNe0+1)
echo foundPNe0 = $foundPNe0
set dimen(statra)=$($(2*$rang+1)*$(2*$rang+1))
define idum 0
do x = -($rang), $rang {
do y = -($rang), $rang {
set
statra[$idum]=$(image[$rx+$x,$ry+$y])
define idum $($idum+1)
}
}
#handles masked pixels
set statrafinal=statra if (statra>-900)
#concat this image onto the growing big file
set woPNe = woPNe concat statrafinal
#end of making
sure we need it
}
#end of found 0 PNe
}
#end of waiting for quotas to be met
}
PNeR
#need to pick random boxes and save their pixels, too
#loop while quota is not med
while {$foundPNeR < $nPNe} {
#pick random coordinates (good
box, for sure)
pickgoodbox rx ry rRA rDEC
#save coordinates
set saveRAr = saveRAr concat $rRA
set saveDECr = saveDECr concat
$rDEC
#collect pixels
set
dimen(statra)=$($(2*$rang+1)*$(2*$rang+1))
define idum 0
do x = $($rx-$rang), $($rx+$rang)
{
do y =
$($ry-$rang), $($ry+$rang) {
set statra[$idum]=$(image[$x,$y])
define idum $($idum+1)
}
}
#handles masked pixels
set statrafinal=statra if
(statra>-800)
#concat this image onto the
growing big vector
set rand = rand concat statrafinal
#iterate counter
define foundPNeR $($foundPNeR+1)
echo foundPNeR = $foundPNeR
#end of foundPNeR < nPNe, quota is not met, loop
}
PNeG
#first, write out files.
!rm $wfound
do w = 0,
$(dimen(saveRAw)-1) { write $wfound "$!(saveRAw[$!w])
$!(saveDECw[$!w]) $!rang" }
!rm $zero
do w = 0, $(dimen(saveRAz)-1) { write $zero
"$!(saveRAz[$!w]) $!(saveDECz[$!w]) $!rang" }
!rm $r
do wr = 0, $(dimen(saveRAr)-1) { write $r
"$!(saveRAr[$!wr]) $!(saveDECr[$!wr]) $!rang" }
#continue with graphing like normal
#graph'n'time(tm)
#histogramming
set wPNeH = histogram(wPNe:bins)
set woPNeH = histogram(woPNe:bins)
set randH = histogram(rand:bins)
#eliminate spike at far edge that only comes from
tails
set wPNeH[$(DIMEN(wPNeH)-1)] =
wPNeH[$(DIMEN(wPNeH)-2)]
set woPNeH[$(DIMEN(woPNeH)-1)] =
woPNeH[$(DIMEN(woPNeH)-2)]
set randH[$(DIMEN(randH)-1)] =
randH[$(DIMEN(randH)-2)]
#and eliminate spike at beginning, too, for
consistency
set wPNeH[0] = wPNeH[1]
set woPNeH[0] = woPNeH[0]
set randH[0] = randH[1]
vecminmax wPNeH wMin wMax
vecminmax randH randMin randMax
if ($wMax > $randMax) {
define lmax $($($wMax+5)*1.1)
} else {
define lmax $($($randMax+5)*1.1)
}
define wM 0
define woM 0
define randM 0
define SIQR 0
# -do med calcs
stats_med wPNe wM SIQR
stats_med woPNe woM SIQR
stats_med rand randM SIQR
#write out median data to file
write "$!v.$!itn.Im$!imag.R$!rang.med"
"$!imag $!rang
$!wM $!woM $!randM "
#modification to graph to .ps and onscreen both
#pretty much don't need graphs, only concerned with
median values (above)
echo graphing...
do gopher = 1,2 {
if ($gopher == 2) {
device
postencap "$!v.$!itn.Im$!imag.R$!rang.ps"
# device x11
echo making
.ps :
}
if ($gopher == 1) {
device
postencap "$!v.$!itn.Im$!imag.R$!rang.ps"
echo
"$!v.$!itn.Im$!imag.R$!rang.ps"
}
era
lim -5 100 0 $lmax
ctype black
box
ctype red
connect bins wPNeH
relocate 55 $($lmax/1.05)
label magenta = 1+ PNe
relocate 55 $($lmax/2)
label med=$wM
ctype green
connect bins woPNeH
relocate 55 $($lmax/1.05 -
$lmax/30)
label green = 0 PNe
relocate 55 $($lmax/2 - $lmax/30)
label med=$woM
ctype blue
connect bins randH
relocate 55 $($lmax/1.05 -
$lmax/15)
label blue = random boxes
relocate 55 $($($lmax/2)-$lmax/15)
label med=$randM
ctype black
xlabel Counts in a pixel
ylabel N occurances
toplabel Image $imag
range $rang $nPNe boxes
}
echo done - $imag $rang:
"$!v.$!itn.Im$!imag.R$!rang.ps"