~/Virgo/FeldPNe/PNe8alpha3/PNe8alpha3.sm
PNe8alpha3
#this program takes our images and makes a
'theoretical' PNe catalog from them
#here's how:
# look at each pixel's flux
# convert flux_v to L_v (d = 16
Mpc)
# convert L_v to L_bol (BC = -0.8)
# multiply L_bol by alpha_2.5 to get N
# set Number of PNe in that field equal to #
where
# if N <= 1 # = N
# if N > 1 # = int(N) + M
#
where M
is:
i.e. decimal portion
#
1 if random number < (int(N) - N)
#
0 if random number > (int(N) - N)
#set up constants we're assuming
#distance, in Mpc
define d 16
#pixel scale, arcseconds^2/pixel
define ps 1.45
#pixel scale, degrees / pixel (from fits
header, in this case)
define dp 4.0277777e-4
#bolometric correction
define BC -0.8
#zero point - surface brightness/count in
mag/arcsec^2
define sbc 29.2
#absolute magnitude of sun (bolometric)
define Mbsun 4.7
#calculate more conversions/constants
#scale/distance: pc^2 / arcmin^2
#using d^2 = (alpha*D/206265)^2
define pc2am2
$($($(60*$d*1000000)/(206265))**2)
#find bolometric apparent mag of sun at 206265
pc
define mbSun $($(5*$(lg(206265))-5 + $Mbsun))
#use this in surface brightness eq. to find
L(bol)/pc^2 per count (w/ zero point)
define Lpc2ct $(10**$($($sbc-$mbSun)/$(-2.5)))
#now find L_sun(bol)/arcmin^2 by combining
Lpc2ct and pc2am2
define Lsunarcmin2 $($Lpc2ct*$pc2am2)
#convert degrees/pixel into square degrees/
(square) pixel
define d2p2 $($dp**2)
#convert deg^2/px(^2) to arcmin^2/px(^2)
define am2p2 $($d2p2*60*60)
#use arcmin^2/px(^2) and L_sun(bol)/arcmin^2
to get L_sun(bol) / px(^2)
define Lsunpx2 $($am2p2*$Lsunarcmin2)
#apply bolometric correction
(m1-m2=-2.5log(f1/f2))
define Lsunpx2BC
$($Lsunpx2*$(10**$($BC/$(-2.5))))
#THIS ^^^^^^^^ is the conversion used
later, it means that each count
# in each pixel is $Lsunpx2BC Solar
Luminosities (bol)
#set up images
#set images = { 3mm.fits 4mm.fits 7mm.fits
Submm.fits }
set images = { Submm.fits }
#set alpha paramters
#set alpha25s = { 23e-9 33e-9 11e-9 }
set alpha25s = { 11e-9 }
#start up loops for images/alphas
foreach imag images {
define file_type fits
image $imag
#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
foreach alpha25 alpha25s {
#time to do
the real work!
PNe8alphaw
}
}
PNe8alphaw
#here's the real work!
#set running total of (ALL) PNe
define currtot 0
#prep vectors to store x,y RA,DEC of PNe
set dimen(PNeRA) = 0
set dimen(PNeDEC) = 0
set dimen(PNex) = 0
set dimen(PNey) = 0
set dimen(PNeRA2) = 0
set dimen(PNeDEC2) = 0
set dimen(PNex2) = 0
set dimen(PNey2) = 0
#start by looping over the entire image, pixel by
pixel
do x = 0, $($NAXIS1-1) {
do y = 0, $($NAXIS2-2) {
#convert with
alpha and L_sun/pixel^2, into "N" PNe
define N
$($(image[$x,$y])*$Lsunpx2BC*$alpha25)
# for
conversion details for ^^^^^^, see top of program
if ($N >
0.5) { echo $N }
define decim
$($N - $(int($N)))
if ($decim
>= $(random(1))) {
define N $($(int($N)+1))
} else {
define N $(int($N))
}
do w = 1, $N {
set PNex = PNex concat $x
set PNey = PNey concat $y
}
if ($N > 1)
{
echo We've got a DOUBLE (or more)!
echo at $x $y we had $N PNe
#write to special vector/file
set PNex2 = PNex2 concat $x
set PNey2 = PNey2 concat $y
}
define currtot
$($currtot + $N)
#end y loop
}
if ($($($x+1)/15) ==
$(INT($($x+1)/15))) {
echo $($x+1) /
$NAXIS1 currtot is $currtot
}
#end x loop
}
#convert PNex,PNey into PNeRA,PNeDEC
set PNeDEC = ((PNey - $CRPIX2 + 1)*$(ABS($CD2_2)) +
$CRVAL2)
set PNeRA = ((-(PNex - $CRPIX1 + 1
))*(ABS($CD1_1)/cosd(PNeDEC)) + $CRVAL1)
set PNeDEC2 = ((PNey2 - $CRPIX2 + 1)*$(ABS($CD2_2))
+ $CRVAL2)
set PNeRA2 = ((-(PNex2 - $CRPIX1 + 1
))*(ABS($CD1_1)/cosd(PNeDEC2)) + $CRVAL1)
#calculate max PNe/pixel, if applicable
define mpp 1
#check for 2+ PNe on same pixel - in PNeRA2 PNeDEC2
list
if ($(dimen(PNeRA2)) > 0 ) {
do wr = 0, $(dimen(PNeRA2)-1) {
write
"8.PNe.$!imag.alpha25$!wrtalpha.RADEC2.dat" $(PNeRA2[$wr])
$(PNeDEC2[$wr])
define mpp
$($mpp + 1)
}
}
define wrtalpha $(int($alpha25 / $(1e-9)))
do wr = 0, $(dimen(PNeRA)-1) {
write
"8.PNe.$!imag.alpha25$!wrtalpha.RADEC.dat" $(PNeRA[$wr]) $(PNeDEC[$wr])
}
#convert PNeRA,PNeDEC to *.reg file for ds9 to
display
set PNeRAs = STRING(PNeRA)
set PNeDECs = STRING(PNeDEC)
set COORD = '"+circle("' + PNeRAs + '","' + PNeDECs
+ '",5) # green"'
#write to file
!rm "8.PNe.$!imag.alpha25.$!wrtalpha.reg"
write "8.PNe.$!imag.alpha25.$!wrtalpha.reg" "#
filename: $!imag"
write "8.PNe.$!imag.alpha25.$!wrtalpha.reg" "#
format: degrees (fk5)"
do wr = 0, $(dimen(PNeRA)-1) {
write
"8.PNe.$!imag.alpha25.$!wrtalpha.reg" $!(COORD[$!wr])
}
write "8.PNe.$!imag.alpha25.$!wrtalpha.reg" "###
Maximum PNe/pixel, anywhere: $!mpp"
#set ds9 image size and zoom for jpg taking
if ('"$!imag"' == '"3mm.fits"') {
define geo 669x880
define zoo 1.00
}
if ('"$!imag"' == '"4mm.fits"') {
define geo 588x798
define zoo 1.00
}
if ('"$!imag"' == '"7mm.fits"') {
define geo 750x950
define zoo 0.500
}
if ('"$!imag"' == '"Submm.fits"') {
define geo 700x900
define zoo 0.500
}
#generate .jpg file of image with overlaid PNe
positions
!rm 8.PNe.$imag.alpha25$wrtalpha.jpg
!"ds9 -fits $!imag -geometry $!geo -zoom $!zoo
-scale limits -50 200 -cmap value 9 0.3 -wcs sky fk5 -wcs skyformat
degrees -region 8.PNe.$!imag.alpha25.$!wrtalpha.reg -saveimage jpeg
8.PNe.$!imag.alpha25$!wrtalpha.jpg -exit"