~/Virgo/FeldPNe/PNe8alpha3/subs.sm


subs
    #subsampling PNe catalogs for alpha_{2.5} = 23e-9
    data 8.PNe.3mm.fits.alpha2523.RADEC.dat
    read { RA3 1 DEC3 2 }
    
    data 8.PNe.4mm.fits.alpha2523.RADEC.dat
    read { RA4 1 DEC4 2 }
    
    data 8.PNe.7mm.fits.alpha2523.RADEC.dat
    read { RA7 1 DEC7 2 }
    
    data 8.PNe.Submm.fits.alpha2523.RADEC.dat
    read { RASub 1 DECSub 2 }
    
    #set contamination fractions to use
    set cs = 0, 10
#    set cs = { 5 }
    
    foreach cr cs {
        
        define c $($cr/10)
        
        #set limits of subsample
        define num3 75
        define num4 11
        define num7 141
        define numSub 34
        
        define w 0
        set dimen(PNeRA) = $($num3+$num4+$num7+$numSub)
        set dimen(PNeDEC) = dimen(PNeRA)
        
        #set up real/contam index for all selected PNe
        set dimen(rea) = dimen(PNeRA)
        #if rea[i] == 1, then i-th value is real, if rea[i] == 0, then i-th value is contam
        
        #do the subsampling
        #for 3
        set select3 = random($(dimen(RA3)))
        sort {select3 RA3 DEC3}
        do i = 0, $($num3-1) {
            set PNeRA[$w] = $(RA3[$i])
            set PNeDEC[$w] = $(DEC3[$i])
            set rea[$w] = 1
            define w $($w+1)
        }
        #DO the contamination fraction for 3
        #read in 3.fits to get header info (i.e. range to choose
        #   random contaminants from)
        define file_type fits
        image 3.fits
        #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
        
        #pick the right number of random coords within this image's
        #   boundaries - in pixels (x,y)
        set randx = ($NAXIS1*(RANDOM(int($($c*$num3)))))
        set randy = ($NAXIS2*(RANDOM(int($($c*$num3)))))
        
        #convert these to RA/DEC
                set randDEC = ((randy - $CRPIX2 + 1)*$(ABS($CD2_2)) + $CRVAL2)
                set randRA = ((-(randx - $CRPIX1 + 1 ))*(ABS($CD1_1)/cosd(randDEC)) + $CRVAL1)
        
        #substitute these back into PNeRA/PNeDEC for appropriate spots
        do i = 0, $(int($c*$num3)-1) {
#            echo substituting in $(randRA[$i]) for $(PNeRA[$($w-$i-1)]), which is the $i th value from the end so far
            set PNeRA[$($w-$i-1)] = $(randRA[$i])
#            echo substituting in $(randDEC[$i]) for $(PNeDEC[$($w-$i-1)]), which is the $i th value from the end so far
            set PNeDEC[$($w-$i-1)] = $(randDEC[$i])
            set rea[$($w-$i-1)] = 0
        }
        
        #for 4
        set select4 = random($(dimen(RA4)))
        sort {select4 RA4 DEC4}
        do i = 0, $($num4-1) {
            set PNeRA[$w] = $(RA4[$i])
            set PNeDEC[$w] = $(DEC4[$i])
            set rea[$w] = 1
            define w $($w+1)
        }
        #DO the contamination fraction for 4
        #read in 4.fits to get header info (i.e. range to choose
        #   random contaminants from)
        define file_type fits
        image 4.fits
        #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
        
#        echo $($c*$num4)
#        echo $(int($($c*$num4)))
        
        #pick the right number of random coords within this image's
        #   boundaries - in pixels (x,y)
        set randx = ($NAXIS1*(RANDOM(int($($c*$num4)))))
        set randy = ($NAXIS2*(RANDOM(int($($c*$num4)))))
        
        #convert these to RA/DEC
                set randDEC = ((randy - $CRPIX2 + 1)*$(ABS($CD2_2)) + $CRVAL2)
                set randRA = ((-(randx - $CRPIX1 + 1 ))*(ABS($CD1_1)/cosd(randDEC)) + $CRVAL1)
        
        #substitute these back into PNeRA/PNeDEC for appropriate spots
        do i = 0, $(int($c*$num4)-1) {
#            echo substituting in $(randRA[$i]) for $(PNeRA[$($w-$i-1)]), which is the $i th value from the end so far
            set PNeRA[$($w-$i-1)] = $(randRA[$i])
#            echo substituting in $(randDEC[$i]) for $(PNeDEC[$($w-$i-1)]), which is the $i th value from the end so far
            set PNeDEC[$($w-$i-1)] = $(randDEC[$i])
            set rea[$($w-$i-1)] = 0
        }
        
        #for 7
        set select7 = random($(dimen(RA7)))
        sort {select7 RA7 DEC7}
        do i = 0, $($num7-1) {
            set PNeRA[$w] = $(RA7[$i])
            set PNeDEC[$w] = $(DEC7[$i])
            set rea[$w] = 1
            define w $($w+1)
        }
        #DO the contamination fraction for 7
        #read in 7.fits to get header info (i.e. range to choose
        #   random contaminants from)
        define file_type fits
        image 7.fits
        #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
        
        #pick the right number of random coords within this image's
        #   boundaries - in pixels (x,y)
        set randx = ($NAXIS1*(RANDOM(int($($c*$num7)))))
        set randy = ($NAXIS2*(RANDOM(int($($c*$num7)))))
        
        #convert these to RA/DEC
                set randDEC = ((randy - $CRPIX2 + 1)*$(ABS($CD2_2)) + $CRVAL2)
                set randRA = ((-(randx - $CRPIX1 + 1 ))*(ABS($CD1_1)/cosd(randDEC)) + $CRVAL1)
        
        #substitute these back into PNeRA/PNeDEC for appropriate spots
        do i = 0, $(int($c*$num7)-1) {
#            echo substituting in $(randRA[$i]) for $(PNeRA[$($w-$i-1)]), which is the $i th value from the end so far
            set PNeRA[$($w-$i-1)] = $(randRA[$i])
#            echo substituting in $(randDEC[$i]) for $(PNeDEC[$($w-$i-1)]), which is the $i th value from the end so far
            set PNeDEC[$($w-$i-1)] = $(randDEC[$i])
            set rea[$($w-$i-1)] = 0
        }
        
        #for Sub
        set selectSub = random($(dimen(RASub)))
        sort {selectSub RASub DECSub}
        do i = 0, $($numSub-1) {
            set PNeRA[$w] = $(RASub[$i])
            set PNeDEC[$w] = $(DECSub[$i])
            set rea[$w] = 1
            define w $($w+1)
        }
        #DO the contamination fraction for Sub
        #read in Sub.fits to get header info (i.e. range to choose
        #   random contaminants from)
        define file_type fits
        image Sub.fits
        #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
        
        #pick the right number of random coords within this image's
        #   boundaries - in pixels (x,y)
        set randx = ($NAXIS1*(RANDOM(int($($c*$numSub)))))
        set randy = ($NAXIS2*(RANDOM(int($($c*$numSub)))))
        
        #convert these to RA/DEC
                set randDEC = ((randy - $CRPIX2 + 1)*$(ABS($CD2_2)) + $CRVAL2)
                set randRA = ((-(randx - $CRPIX1 + 1 ))*(ABS($CD1_1)/cosd(randDEC)) + $CRVAL1)
        
        #substitute these back into PNeRA/PNeDEC for appropriate spots
        do i = 0, $(int($c*$numSub)-1) {
#            echo substituting in $(randRA[$i]) for $(PNeRA[$($w-$i-1)]), which is the $i th value from the end so far
            set PNeRA[$($w-$i-1)] = $(randRA[$i])
#            echo substituting in $(randDEC[$i]) for $(PNeDEC[$($w-$i-1)]), which is the $i th value from the end so far
            set PNeDEC[$($w-$i-1)] = $(randDEC[$i])
            set rea[$($w-$i-1)] = 0
        }
        
        #echo $c
        #echo PNeFAKE.C$cr.dat
        !rm PNeFAKE.C$cr.dat
        write PNeFAKE RA DEC real
        do i = 0, $(dimen(PNeRA)-1) {
            write PNeFAKE.C$cr.dat $(PNeRA[$i]) $(PNeDEC[$i]) $(rea[$i])
        }
    }