/Virgo/FeldPNe/PNe8/PNe.sm
OR
/Virgo/FeldPNe/vat/PNe8/PNe.sm

PNe
    #latest version 8! only three histograms: With PNe, Without PNe, random
    #   should eliminate some bias
    #
    #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 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)
   
    #define images/ranges
    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 = { 3m.fits }
    #set raO = { 10 15 20 }
   
    data PNe.param
    read { SZ 1 SO 2 SR 3 }
    define SpaceZero $(SZ[0])
    define SpaceOne $(SO[0])
    #define SpaceTwo $(ST[0])
    define SpaceRand $(SR[0])
   
    #leave this zero
    define num 0
   
    #big images/ranges loops
    foreach imag imO {
       
        define file_type fits
        image $imag
       
        foreach rang raO {
           
            #run these
            #prep data
            PNew
           
            #set number of boxes here
            #define reals 100
            #for now it's based on number of PNe in this field
            define reals 0
            #establishing number of PNe in this field:
            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 && $rx > 0 && $ry < $NAXIS2 && $ry > 0) {
                    define reals $($reals + 1)
                }
            }
           
            #run boxes with PNe
            PNe1
           
            #run boxes without PNe
            PNe0
           
            #run random boxes
            PNer
           
            #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

   
    #read in PNe datafile
    data PNe.dat
    read { field 1 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(wPNe) = 0
    set DIMEN(woPNe) = 0
    set DIMEN(rand) = 0


PNe1
    #random boxes that contain one (or more) PN
   
    #file to save box coordinates in and compare to later
    define wfound "$!v.$!itn.Im$!imag.R$!rang.wfound.dat"
   
        set dimen(saveRAw) = $reals
        set dimen(saveDECw) = $reals
   
    #loop for as many boxes as we want
    define makesure 0
    define num 0
    while {$num < $reals} {
        #generate random (x,y) coordinates
                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) {
            #found at least one PN in this box
           
            #check for other boxes too close, if SpaceOne is one
            set delRAw = $rRA - saveRAw
                        set delDECw = $rDEC - saveDECw
           
             #check for SpaceOne
                        if ($SpaceOne == 1) {
                                set counterp = ( (abs(delRAw)<=$($rangeRD)) && (abs(delDECw)<=$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
                #save coords and save pixels
                set saveRAw[$num] = $rRA
                                set saveDECw[$num] = $rDEC
                                define num $($num + 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])
                                                #handles masked pixels
                                                set statrafinal=statra if (statra>-900)
                                                define idum $($idum+1)
                                        }
                                }
               
                                #concat this image onto the growing big file
                                set wPNe = wPNe concat statrafinal
               
            #end of checking that other boxes aren't too close   
            } ELSE {
                define makesure $($makesure + 1)
            }
       
        #end of is 1+ PN found
        } ELSE {
            define makesure $($makesure + 1)
        }
       
        #check makesure - don't let it look more than 10000 times #_ boxes
         if ($makesure > $($reals*10000)) {
                        echo There is trouble - Couldn't get enough boxes in $($reals*10000) trials... you suck
                        define numz $num
                        define num $reals
                }
   
    #end while still doing more boxes, while $num < $reals, loop
    }
   
     #write out data file
        !rm $wfound
        do w = 0, $(dimen(saveRAw)-1) { write $wfound "$!(saveRAw[$!w])  $!(saveDECw[$!w]) $!rang" }

        echo 33% - did 1+ PNe/box WITH optional spacing
   
   
PNe0
    #checking for boxes WITHOUT PNe in them. # boxes = $reals
   
    #file to save box coordinates in and compare to later
    define zero "$!v.$!itn.Im$!imag.R$!rang.0found.dat"
   
    set dimen(saveRAz) = $reals
    set dimen(saveDECz) = $reals
   
    define makesure 0
   
    define num 0

    while {$num < $reals} {
        #pick random (x,y)
        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) {
           
            #make sure no other boxes close to here
           
            set delRAz = $rRA - saveRAz
            set delDECz = $rDEC - saveDECz
           
            #check for SpaceZero
            if ($SpaceZero == 1) {
                set counterp = ( (abs(delRAz)<=$($rangeRD)) && (abs(delDECz)<=$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])
                        #handles masked pixels
                           set statrafinal=statra if (statra>-900)
                        define idum $($idum+1)
                    }
                }
               
                #concat this image onto the growing big file
                set woPNe = woPNe 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 There is trouble - Couldn't get enough boxes in $($reals*10000) trials... you suck
            define numz $num
            define num $reals
        }
    #end while $num < $reals loop for overall 0 PNe loop
    }

    #write out data file
    !rm $zero
    do w = 0, $(dimen(saveRAz)-1) { write $zero "$!(saveRAz[$!w])  $!(saveDECz[$!w]) $!rang" }
       
    echo 66% - did 0 PNe/box WITH optional spacing       
   
   
   
PNer
    #random boxes
   
    #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.Rfound.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)
       
        #make sure no other boxes close to here
           
        set delRA = $rRA - saveRAr
        set delDEC = $rDEC - saveDECr
       
        #check for SpaceRand
        if ($SpaceRand == 1) {
            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])
                    #handles masked pixels
                       set statrafinal=statra if (statra>-900)
                    define idum $($idum+1)
                }
           
            }
           
            #concat this image onto the growing big file
            set rand = rand concat statrafinal
        #end of if sum == 0 loop
        } ELSE {
            define makesurer $($makesurer + 1)
        }
       
        if ($makesurer > $($reals*10000)) {
            echo There is trouble - Couldn't get enough boxes in $($reals*10000) trials... you suck
            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)
   
    #histogramming
    set wPNeH = histogram(wPNe:bins)
    set woPNeH = histogram(woPNe:bins)
    set randH = histogram(rand:bins)
   
    #eliminate spike at 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)]
   
   
    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    $reals boxes
    }
    echo done - $imag $rang: "$!v.$!itn.Im$!imag.R$!rang.ps"