/packmule_raid/steven/PNe7/PNe.sm


PNe
    #latest, v.7 will allow user to select which methods to space boxes apart
    #     in the file PNe.param
    #v.7 also reads in i.i to establish itn, iteration, and will do many of
    #    them, to establish error bars
    #v.7 also scrambles PNe catalog order before using them in On(e)PNe and
    #    TwoPNe methods, so as to make the selection of pairs/non-overlapping
    #    individuals totally random
    #
    #this version (7) runs a COMPLETE iteration each time it's called, i.e. on each
    #    node. goes through every imO and raO, then stops
    #
    #latest, v.6 will mandate that 0PNe boxes are a half-box-width apart, too
    #
    #this version implements a new 2+ PNe selection technique and makes sure
    #    they're spaced apart!
    #
    ##generates 4 histograms:
    ## 1) N/I pixels for boxes with 2 or more PNe in them, randomly selected
    ## 2) N/I pixels for boxes around PNe
    ## 3) N/I pixels for boxes WITHOUT PNe in them
    ## 4) N/I pixels for totally randomly selected boxes
    
    #read in PNe datafile
    data PNe.dat
    read { field 1.s num 2 RAh 3 RAm 4 RAs 5 DECd 6 DECam 7 DECas 8 m5007 9 notes 10.s }

    #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)
    
    #set imO = { 3.fits 3m.fits 4.fits 4m.fits 7.fits 7m.fits Sub.fits Subm.fits }
    #set raO = { 10 15 20 25 30 35 40 45 50 55 60 }
    
    #temp for testing
    set imO = { Core.fits Corem.fits FCJ.fits FCJm.fits }
#    set imO = { 4.fits 4.fits }
    set raO = { 1 3 5 7 9 11 13 15 17 20 25 30 35 40 45 50 55 60 }
    
    #read in whether or not to use spacing for boxes
    data PNe.param
    read {SZ 1 SO 2 ST 3 SR 4}
    define SpaceZero $(SZ[0])
    define SpaceOne $(SO[0])
    define SpaceTwo $(ST[0])
    define SpaceRand $(SR[0])
    
    #big loops for images and ranges
    foreach rang raO {
        foreach imag imO {
    
            define file_type fits
            image $imag
            
            #run everything!
            #prep data/vectors
            PNew
            
            #run 1 box on each PNe
            PNe1
            
            #run 2 PNe/box
            PNe2
            
            #run 0 PNe/box (with #1PNe/#2PNe)
            PNe31
            PNe32
            
            #run totally random (with #1PNe/#2PNe)
            PNe41
            PNe42
            
            #run grapher
            PNeg
            
        #end of images loop
        }
    #end of ranges loop
    }    

PNew
    #this does the real work and is called by PNe, which is just loops
    
    #define as fits and read in important header info

    #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
    
    #prep for graphing inline
    set bins = -5, 100, 0.5

#    #done at top now    
#    #read in PNe datafile
#    data PNe.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 DEC = DECd + (DECam/60) + (DECas/3600)
    set RA = 15 * (RAh + (RAm/60) + (RAs/3600))
    define rangeRD (($rang)*ABS($CD1_1))
    #^these conversions are checked and good. duh.^#
    
    #randomise RA/DEC for use in this iteration
    #set scram = RANDOM($(DIMEN(RA)))
    #sort {scram RA DEC}
    
    set DIMEN(onPNe) = 0
    set DIMEN(offPNe1) = 0
    set DIMEN(offPNe2) = 0
    set DIMEN(rand1) = 0
    set DIMEN(rand2) = 0
    set DIMEN(xkeepo) = 1    
    set DIMEN(ykeepo) = 1
    define reals 0
    define num 0
    

PNe1
    #(second histogram, done first)
    #second: place box on each PNe and calculate average light in it
    
    #file to save box coordinates in and compare to later
    define on "$!v.$!itn.Im$!imag.R$!rang.1found.dat"
    !rm $on
    
    define reals 0
    
    do i = 0, $($(DIMEN(RA)-1)), 1 {
        #convert RA/DEC back to pixels
        #(this is right)
        define ry $(INT($CRPIX2 + ($($(DEC[$i])-$CRVAL2)/$(ABS($CD2_2)))))
        define rx $(INT($($CRPIX1 - 1 - $($($(RA[$i] - $CRVAL1)*cosd($(DEC[$i])))/$(ABS($CD1_1))))))
                    
        #check to make sure PNe is inside this image's bounds with room
        #                to take a region around it
        if ($rx < $NAXIS1-$rang && $rx > $rang && $ry < $NAXIS2-$rang && $ry > $rang) {

            #if SpaceOne == 1, DO check spacing, if SpaceOne == 0, do NOT check spacing
            if ($SpaceOne == 1) {
                #check if too close to any other onPNe boxes
                set delxon = $rx - xkeepo
                set delyon = $ry - ykeepo
                set counteron = ((abs(delxon)<=$($rang)) && (abs(delyon)<=$($rang)))?1:0
            } ELSE {
                set counteron = {0}
            }
            if ($(sum(counteron)) == 0) {
                #store these coords
                set xkeepo = xkeepo concat $rx
                set ykeepo = ykeepo concat $ry
                        
                #keep track of how many used for later
                define reals $($reals + 1)
                
                #collecting from a box around PNe 2*rang+1 x 2*rang+1
                    
#                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 onPNe = onPNe concat statrafinal

            #end of if no other box centers in this box
            }
        #end of making sure PNe is inside this image loop
        }
    #end of checking through catalog for PNe loop
    }
    
    #write out stored coords
    #convert to RA/DEC first
    set DECon = ((ykeepo - $CRPIX2 + 1)*$(ABS($CD2_2)) + $CRVAL2)
    set RAon = ((-(xkeepo - $CRPIX1 + 1 ))*(ABS($CD1_1)/cosd(DECon)) + $CRVAL1)
    
    #header first
    write $on "#RA_coord     DEC_coord     range(px)"
    
    do wrtro = 1, $(DIMEN(RAon)-1) {
        write $on "$!(RAon[$!wrtro])  $!(DECon[$!wrtro]) $!rang"
    }
    
    
    echo 25% - did box on each PNe
    

PNe2    
    #'first' part is done second
    #first part: find all possible boxes, eliminate those too close together if SpaceTwo == 1
    set DIMEN(twoPNe) = 0
    define num 0
    set DIMEN(RAkeep) = $($reals+1)
    set DIMEN(DECkeep) = $($reals+1)

    #variable to keep track of how many boxes found with 2+ PNe
    define nb 0
    
    #file to save box coordinates in and compare to later
    define g "$!v.$!itn.Im$!imag.R$!rang.2found.dat"
    !rm $g
    
    #figure out how many distances will be calculated
    define N $(DIMEN(RA))
    define tot 0
    do i = 1,$N {
        define tot $($tot + $i)
    }
    echo distances to be calculated: $tot
    
    #prep for data
    set DIMEN(xbox)=$tot
    set DIMEN(ybox)=$tot    
    
    #start by IDing all PNe separated by < 2*range + 1
    
    #set up 2 loops to get distances between all PNe in catalog table6.dat
    do L1 = 0, $(DIMEN(RA)-1) {
        do L2 = $($L1 + 1), $(DIMEN(RA)-1) {
            
            #get coordinates of PN #L1 and PN #L2 into x & y
            #note, don't have to be INT values
            define yL1 $(($CRPIX2 + ($($(DEC[$L1])-$CRVAL2)/$(ABS($CD2_2)))))
            define xL1 $(($($CRPIX1 - 1 - $($($(RA[$L1] - $CRVAL1)*cosd($(DEC[$L1])))/$(ABS($CD1_1))))))
            define yL2 $(($CRPIX2 + ($($(DEC[$L2])-$CRVAL2)/$(ABS($CD2_2)))))
            define xL2 $(($($CRPIX1 - 1 - $($($(RA[$L2] - $CRVAL1)*cosd($(DEC[$L2])))/$(ABS($CD1_1))))))
            
            #check to see if both PNe L1 and L2 are in this field
            if ($xL1 < $NAXIS1-$rang && $xL1 > $rang && $yL1 < $NAXIS2-$rang && $yL1 > $rang && $xL2 < $NAXIS1-$rang && $xL2 > $rang && $yL2 < $NAXIS2-$rang && $yL2 > $rang) {
                            
                #check if distance between PN #L1 and PN #L2 is < 2*range+1
                if ($(abs($yL1 - $yL2)) < $(2*$rang+1) && $(abs($xL1 - $xL2)) < $(2*$rang+1)) {

                    #if SpaceTwo == 1, DO check spacing, if SpaceTwo == 0, do NOT check spacing
                    if ($SpaceTwo == 1) {
                        #check if any other boxes are too close to this one
                        set delx = $($($xL1 + $xL2)/2) - xbox
                        set dely = $($($yL1 + $yL2)/2) - ybox
                        set counterr = ((abs(delx)<=$($rang)) && (abs(dely)<=$($rang)))?1:0
                    } ELSE {
                        set counterr = {0}
                    }
                    if ($(sum(counterr)) == 0) {

                        #--> found a close PN pair in this image!
                        #:  store coords, iterate nb
                        #make a box halfway between these
                        set xbox[$nb] = $($($xL1 + $xL2)/2)
                        set ybox[$nb] = $($($yL1 + $yL2)/2)
                        
                        #store values of light in this region
                        set dimen(statra)=$($(2*$rang+1)*$(2*$rang+1))
#                        define idum 0
#                        do x = -($rang), $rang {
#                            do y = -($rang), $rang {
#                                set statra[$idum]=$(image[$(xbox[$nb])+$x,$(ybox[$nb])+$y])
#                                define idum $($idum+1)
#                            }
#                        }
#                        
#                        #handles masked pixels
#                           set statrafinal=statra if (statra>-900)
#                        
                        #concat this image to the growing big file
                        set twoPNe = twoPNe concat statrafinal
                        
                        #iterate nb
                        define nb $($nb + 1)
                        
                    #end of if too close to another box loop
                    }
                #end if distance < whatever loop
                }
            #end if both PNe are in this field loop
            }
        #end looping through 2nd PN loop
        }
    #end looping through 1st PN loop
    }
    
    #output these box coords to file
    #header first
    write $g "#RA_coord     DEC_coord     range(px)"
    
    #get xbox and ybox into RAbox and DECbox
    set DECbox = ((ybox - $CRPIX2 + 1)*$(ABS($CD2_2)) + $CRVAL2)
    set RAbox = ((-(xbox - $CRPIX1 + 1 ))*(ABS($CD1_1)/cosd(DECbox)) + $CRVAL1)
    
    #figure out what x,y = 0,0 is in RA/DEC
    define DECzero $($(0 - $CRPIX2 + 1)*$(ABS($CD2_2)) + $CRVAL2)
    define RAzero $($(-(0 - $CRPIX1 + 1 ))*(ABS($CD1_1)/cosd($DECzero)) + $CRVAL1)

    #write out coords that are not x,y = 0,0, but write in RA,DEC
    do wrtr = 0, $($tot-1) {
        if ($(RAbox[$wrtr]) != $RAzero && $(DECbox[$wrtr]) != $DECzero) {
            #do the write
            write $g "$!(RAbox[$!wrtr])  $!(DECbox[$!wrtr]) $!rang"
        }
    }
    
    
    #end of checking through random boxes with 2+ PNe section
    
    
    echo 50% - did 2+ PNe/box WITH 'range'-spaced boxes

PNe31
        #FOR ONE BOX/PNE nUMBERS
    #third part: random boxes that have 0 PNe in them AND (with v.6) are at least a half-box width (rang) away from any other box with 0 PNe
        #    spaced IFF SpaceZero == 1

        #file to save box coordinates in and compare to later
        define h "$!v.$!itn.Im$!imag.R$!rang.0found1.dat"

        set dimen(saveRAz) = $reals
        set dimen(saveDECz) = $reals

        define makesure 0

        define num 0
        while {$num < $reals} {
                #generate random coordinates as before
                define rx $(($(RANDOM(1))*($NAXIS1-(2*$rang)))+$rang)
                define ry $(($(RANDOM(1))*($NAXIS2-(2*$rang)))+$rang)

                #convert to RA/DEC
                define rDEC  (($ry - $CRPIX2 + 1)*$(ABS($CD2_2)) + $CRVAL2)
                define rRA $($(-($rx - $CRPIX1 + 1 ))*$(ABS($CD1_1)/cosd($rDEC)) + $CRVAL1)

                #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) {

                        #check for SpaceZero
                        if ($SpaceZero == 1) {
                
                            #make sure no other boxes close to here
                
                            set delRA = $rRA - saveRAz
                            set delDEC = $rDEC - saveDECz
                                set counterp = ( (abs(delRA)<=$($rangeRD)) && (abs(delDEC)<=$rangeRD) )? 1 : 0
                        } ELSE {
                                set counterp = {0}
                        }
                        if ($(sum(counterp)) == 0) {

                                #now we have a box that contains 0 and is far enough from other boxes
                                set saveRAz[$num] = $rRA
                                set saveDECz[$num] = $rDEC
                                define num $($num + 1)

                                #collecting from a box around PNe 2*rang+1 x 2*rang+1

                                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 offPNe1 = offPNe1 concat statrafinal
                
                        #end of making sure no other boxes too close
                        } ELSE {
                                define makesure $($makesure + 1)
                        }

                #end of if no PNe found
                } ELSE {
                        define makesure $($makesure + 1)
                }

                #check to make sure it doesn't loop too long
                if ($makesure > $($reals*1000)) {
                        echo Range = $rang is too big to get 0 PNe/box without overlap - histogram will not be normalised
                        define numz $num
                        define num $reals
                }
        #end while $num < $reals loop for overall 0 PNe loop
        }

        #write out data file
        !rm $h
        do w = 0, $(dimen(saveRAz)-1) { write $h "$!(saveRAz[$!w])  $!(saveDECz[$!w]) $!rang" }

        echo 75% - did 0 PNe/box WITH optional spacing

    

PNe32
    #FOR 2+PNE/BOX NUMBERS
    #third part: random boxes that have 0 PNe in them AND (with v.6) are at least a half-box width (rang) away from any other box with 0 PNe
    #    spaced IFF SpaceZero == 1

    #file to save box coordinates in and compare to later
    define h "$!v.$!itn.Im$!imag.R$!rang.0found2.dat"
    
    set dimen(saveRAz) = $nb
    set dimen(saveDECz) = $nb
    
    define makesure 0
    
    define num 0
    while {$num < $nb} {
        #generate random coordinates as before
        define rx $(($(RANDOM(1))*($NAXIS1-(2*$rang)))+$rang)
        define ry $(($(RANDOM(1))*($NAXIS2-(2*$rang)))+$rang)
        
        #convert to RA/DEC
        define rDEC  (($ry - $CRPIX2 + 1)*$(ABS($CD2_2)) + $CRVAL2)
        define rRA $($(-($rx - $CRPIX1 + 1 ))*$(ABS($CD1_1)/cosd($rDEC)) + $CRVAL1)
        
        #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) {
            
            #check for SpaceZero
            if ($SpaceZero == 1) {
                #make sure no other boxes close to here
                
                set delRA = $rRA - saveRAz
                set delDEC = $rDEC - saveDECz
                
                set counterp = ( (abs(delRA)<=$($rangeRD)) && (abs(delDEC)<=$rangeRD) )? 1 : 0
            } ELSE {
                set counterp = {0}
            }
            if ($(sum(counterp)) == 0) {
                
                #now we have a box that contains 0 and is far enough from other boxes
                set saveRAz[$num] = $rRA
                set saveDECz[$num] = $rDEC
                define num $($num + 1)
                
                #collecting from a box around PNe 2*rang+1 x 2*rang+1
                    
                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 offPNe2 = offPNe2 concat statrafinal
                
            #end of making sure no other boxes too close
            } ELSE {
                define makesure $($makesure + 1)
            }

        #end of if no PNe found
        } ELSE {
            define makesure $($makesure + 1)
        }
        
        #check to make sure it doesn't loop too long
        if ($makesure > $($reals*1000)) {
            echo Range = $rang is too big to get 0 PNe/box without overlap - histogram will not be normalised
            define numz $num
            define num $reals
        }
    #end while $num < $reals loop for overall 0 PNe loop
    }

    #write out data file
    !rm $h
    do w = 0, $(dimen(saveRAz)-1) { write $h "$!(saveRAz[$!w])  $!(saveDECz[$!w]) $!rang" }
        
    echo 75% - did 0 PNe/box WITH optional spacing

PNe41    
    #fourth part: random boxes, spaced IFF SpaceRand == 1
    
    #initialise storage vars/vects
    set dimen(saveRAr) = $reals
    set dimen(saveDECr) = $reals
    
    #file to save box coordinates in and compare to later
    define r "$!v.$!itn.Im$!imag.R$!rang.Rfound1.dat"
    
    define numr 0
    define makesurer 0
    
    while {$numr < $reals} {
        
        #generate random coordinates as before
        define rx $(($(RANDOM(1))*($NAXIS1-(2*$rang)))+$rang)
        define ry $(($(RANDOM(1))*($NAXIS2-(2*$rang)))+$rang)
        
        #convert to RA/DEC
        define rDEC  (($ry - $CRPIX2 + 1)*$(ABS($CD2_2)) + $CRVAL2)
        define rRA $($(-($rx - $CRPIX1 + 1 ))*$(ABS($CD1_1)/cosd($rDEC)) + $CRVAL1)
        
        #check for SpaceRand
        if ($SpaceRand == 1) {
            #make sure no other boxes close to here
            
            set delRA = $rRA - saveRAr
            set delDEC = $rDEC - saveDECr
            set counterp = ( (abs(delRA)<=$($rangeRD)) && (abs(delDEC)<=$rangeRD) )? 1 : 0
        } ELSE {
            set counterp = {0}
        }
        if ($(sum(counterp)) == 0) {
            
            set saveRAr[$numr] = $rRA
            set saveDECr[$numr] = $rDEC
            define numr $($numr + 1)
            
            #collecting from a box around PNe 2*rang+1 x 2*rang+1
        
            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 rand1 = rand1 concat statrafinal
        #end of if sum == 0 loop
        } ELSE {
            define makesurer $($makesurer + 1)
        }
        
        if ($makesurer > $($reals*1000)) {
            echo Range = $rang is too big to get Random boxes without overlap - histogram will not be normalised
            define numrr $numr
            define numr $reals
        }
    #end of checking through random boxes
    }
    
    #write data file out
    !rm $r
    do wr = 0, $(dimen(saveRAr)-1) { write $r "$!(saveRAr[$!wr])  $!(saveDECr[$!wr]) $!rang" }

    
    echo 95% - did random boxes

PNe42
    #fourth part: random boxes, spaced IFF SpaceRand == 1
    
    #initialise storage vars/vects
    set dimen(saveRAr) = $nb
    set dimen(saveDECr) = $nb
    
    #file to save box coordinates in and compare to later
    define r "$!v.$!itn.Im$!imag.R$!rang.Rfound2.dat"
    
    define numr 0
    define makesurer 0
    
    while {$numr < $nb} {
        
        #generate random coordinates as before
        define rx $(($(RANDOM(1))*($NAXIS1-(2*$rang)))+$rang)
        define ry $(($(RANDOM(1))*($NAXIS2-(2*$rang)))+$rang)
        
        #convert to RA/DEC
        define rDEC  (($ry - $CRPIX2 + 1)*$(ABS($CD2_2)) + $CRVAL2)
        define rRA $($(-($rx - $CRPIX1 + 1 ))*$(ABS($CD1_1)/cosd($rDEC)) + $CRVAL1)
        
        #check for SpaceRand
        if ($SpaceRand == 1) {
            #make sure no other boxes close to here
                
            set delRA = $rRA - saveRAr
            set delDEC = $rDEC - saveDECr
            
            set counterp = ( (abs(delRA)<=$($rangeRD)) && (abs(delDEC)<=$rangeRD) )? 1 : 0
        } ELSE {
            set counterp = {0}
        }
        if ($(sum(counterp)) == 0) {
            
            set saveRAr[$numr] = $rRA
            set saveDECr[$numr] = $rDEC
            define numr $($numr + 1)
            
            #collecting from a box around PNe 2*rang+1 x 2*rang+1
        
            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 rand2 = rand2 concat statrafinal
        #end of if sum == 0 loop
        } ELSE {
            define makesurer $($makesurer + 1)
        }
        
        if ($makesurer > $($reals*1000)) {
            echo Range = $rang is too big to get Random boxes without overlap - histogram will not be normalised
            define numrr $numr
            define numr $reals
        }
    #end of checking through random boxes
    }
    
    #write data file out
    !rm $r
    do wr = 0, $(dimen(saveRAr)-1) { write $r "$!(saveRAr[$!wr])  $!(saveDECr[$!wr]) $!rang" }

    
    echo 95% - did random boxes
    echo running PNeg to graph them
    
    

PNeg
    #graph'n'time(tm)
    echo starting graphing
    
    #histogramming
    set onPNeH = histogram(onPNe:bins)
    set offPNeH1 = histogram(offPNe1:bins)
    set offPNeH2 = histogram(offPNe2:bins)
    set randH1 = histogram(rand1:bins)
    set randH2 = histogram(rand2:bins)
    set twoPNeH = histogram(twoPNe:bins)
    
    #eliminate spike at edge that only comes from tails
    set onPNeH[$(DIMEN(onPNeH)-1)] = onPNeH[$(DIMEN(onPNeH)-2)]
    set offPNeH1[$(DIMEN(offPNeH1)-1)] = offPNeH1[$(DIMEN(offPNeH1)-2)]
    set offPNeH2[$(DIMEN(offPNeH2)-1)] = offPNeH2[$(DIMEN(offPNeH2)-2)]
    set randH1[$(DIMEN(randH1)-1)] = randH1[$(DIMEN(randH1)-2)]
    set randH2[$(DIMEN(randH2)-1)] = randH2[$(DIMEN(randH2)-2)]
    set twoPNeH[$(DIMEN(twoPNeH)-1)] = twoPNeH[$(DIMEN(twoPNeH)-2)]
    
    
    vecminmax twoPNeH twoMin twoMax
    vecminmax randH1 randMin randMax
    if ($twoMax > $randMax) {
        define lmax $($($twoMax+5)*1.1)
        } else {
        define lmax $($($randMax+5)*1.1)
        }
    
    define onM 0
    define offM1 0
    define offM2 0
    define randM1 0
    define randM2 0
    define twoM 0
    define SIQR 0
    
    # -do med calcs
    stats_med onPNe onM SIQR
    stats_med offPNe1 offM1 SIQR
    stats_med offPNe2 offM2 SIQR
    stats_med rand1 randM1 SIQR
    stats_med rand2 randM2 SIQR
    stats_med twoPNe twoM SIQR
    
    #write out median data to file
    write "$!v.$!itn.Im$!imag.R$!rang.med" "$!imag    $!rang    $!twoM    $!onM    $!offM1    $!offM2    $!randM1    $!randM2 "
    #modification to graph to .ps and onscreen both
    
    #pretty much don't need graphs, only concerned with median values (above)
    
    do gopher = 1,2 {
        if ($gopher == 2) {
            device postencap "$!v.$!itn.Im$!imag.R$!rang.ps"
            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 magenta
        connect bins twoPNeH
        relocate 55 $($lmax/1.05)
        label magenta = 2+ PNe
        relocate 55 $($lmax/2)
        label med=$twoM
        
        ctype red
        connect bins onPNeH
        relocate 55 $($lmax/1.05 - $lmax/30)
        label red = around PNe
        relocate 55 $($lmax/2 - $lmax/30)
        label med=$onM
        
        ctype green
        connect bins randH1
        relocate 55 $($lmax/1.05 - $lmax/15)
        label green = random boxes
        relocate 55 $($($lmax/2)-$lmax/15)
        label med=$randM1
        
        ctype blue
        connect bins offPNeH1
        relocate 55 $($lmax/1.05 - $lmax/10)
        label blue = away from PNe
        relocate 55 $($($lmax/2)-$lmax/10)
        label med=$offM1
        
        ctype black
        xlabel Counts in a pixel
        ylabel N occurances
        toplabel Image $imag   range $rang  
    }
    echo done - $imag $rang: "$!v.$!itn.Im$!imag.R$!rang.ps"