/Virgo/FeldPNe/PNe8alpha/PNe8alpha.sm
OR
/Virgo/FeldPNe/vat/PNe8alpha/PNe8alpha.sm

PNe8alpha
    #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 images
    set images = { 3mm.fits 4mm.fits 7mm.fits Submm.fits }
        
    #set alpha paramters
    set alpha25s = { 23e-9 33e-9 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])*$(8.527e2)*$alpha25*2.09)
            #                                                                             ^BC^
            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"