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

PNeError
    #two part program
    #   -first to compute the values for the values for the error bars
    #   -second to graph
    #      -the average values for each image/range over all iterations
    #      -the error bars!!
   
    #first set up vectors of all possibilities
    set im = { 3.fits 3m.fits 4.fits 4m.fits 7.fits 7m.fits Sub.fits Subm.fits }
    data v.v
    read { vnum 1 }
    define v $(vnum[0])
    data i.i
    read { inum 1 }
    define itt $(inum[0])
    set it = 1,$itt
    set ra = 10,60,5
   
    #temp for testing
    set im = {7.fits}
    set ra = 10,50,5
   
    #separate main runings, first by imagee
    foreach imag im {
       
        #prep vectors for graphing
        set dimen(wGm) = 0
        set dimen(wGs) = 0
        set dimen(wGu) = 0
        set dimen(wGl) = 0
        set dimen(woGm) = 0
        set dimen(woGs) = 0
        set dimen(woGu) = 0
        set dimen(woGl) = 0
        set dimen(randGm) = 0
        set dimen(randGs) = 0
        set dimen(randGu) = 0
        set dimen(randGl) = 0
        set dimen(wMt) = 0
        set dimen(woMt) = 0
        set dimen(randMt) = 0
       
        #next loop is by image
        foreach rang ra {
       
            #call computational program - it loops through ranges
            PNeErrorc
           
        #end of ranges loop
        }
       
        #call graphing program - now all ranges for imag have been done
        PNeErrorg
       
        #don't ask why but it's happier if you run it twice
        #as time permits, leave the next line un-commented-out
        PNeErrorg
       
    #end of images loop
    }
    echo converting to gifs/uploading to the web and Virgo/Web folder...
    qc
    echo done

PNeErrorc
    #this program computes error bars for the PNe7 run, based on .med files
    #  which should be labeled 7.#.Im3.fits.R10.med
    #    where # is iteration, 1 to ____
    #    3.fits file name, 3,3m,4,4m,7,7m,Sub,Subm
    #    R10 is range, 10 to 60, by 5s
    #  file will look like:
    #    write "$!v.$!itn.Im$!imag.R$!rang.med" "$!imag    $!rang    $!wM    $!woM    $!randM "
   
           
            set DIMEN(wMt) = 0
            set DIMEN(woMt) = 0
            set DIMEN(randMt) = 0
            #innermost loop is by iteration number
            foreach itn it {
                #clear variables
                set dimen(wMv) = 0
                set dimen(woMv) = 0
                set dimen(randMv) = 0
               
                #read in data IF file exists
                if (is_file($v.$itn.Im$imag.R$rang.med)) {
                    data $v.$itn.Im$imag.R$rang.med
                    read { imchk 1.s rachk 2 wMv 3 woMv 4 randMv 5 }
                   
                    #store values in temp vectors
                    set wMt = wMt concat $(wMv[0])
                    set woMt = woMt concat $(woMv[0])
                    set randMt = randMt concat $(randMv[0])
                }
               
            #end of iterations loop
            }
           
            #all iterations for one range/one image - stats time
           
            #finding sigma:
            stats wMt wMean wSigma wK
            stats woMt woMean woSigma woK
            stats randMt randMean randSigma randK
           
            #stats_med does not output quartiles so copy/do manually
            sort {wMt}
            define wMed (wMt[DIMEN(wMt)/2])
            define wLQ (wMt[0.25*DIMEN(wMt)])
            define wUQ (wMt[0.75*DIMEN(wMt)])
            sort {woMt}
            define woMed (woMt[DIMEN(woMt)/2])
            define woLQ (woMt[0.25*DIMEN(woMt)])
            define woUQ (woMt[0.75*DIMEN(woMt)])
            sort {randMt}
            define randMed (randMt[DIMEN(randMt)/2])
            define randLQ (randMt[0.25*DIMEN(randMt)])
            define randUQ (randMt[0.75*DIMEN(randMt)])
           
            #write results to data file
            !rm "$!v.0.Im$!imag.R$!rang.stats"
            #writerand1Gm header
            write "$!v.0.Im$!imag.R$!rang.stats" "#twoPNe onPNe offPNe1 offPNe2 rand1 rand2 - data rows go: (1)Mean, (2)Sigma, (3)Kurtosis, (4)Median, (5)Lower Quartile, (6)Upper Quartile"
            write "$!v.0.Im$!imag.R$!rang.stats" $wMean $woMean $randMean
            write "$!v.0.Im$!imag.R$!rang.stats" $wSigma $woSigma $randSigma
            write "$!v.0.Im$!imag.R$!rang.stats" $wK $woK $randK
            write "$!v.0.Im$!imag.R$!rang.stats" $wMed $woMed $randMed
            write "$!v.0.Im$!imag.R$!rang.stats" $wLQ $woLQ $randLQ
            write "$!v.0.Im$!imag.R$!rang.stats" $wUQ $woUQ $randUQ
           
            #store variables for graphing soon
            set wGm = wGm CONCAT $wMean
            set wGs = wGs CONCAT $wSigma
            set wGu = wGu CONCAT $wUQ
            set wGl = wGl CONCAT $wLQ
            set woGm = woGm CONCAT $woMean
            set woGs = woGs CONCAT $woSigma
            set woGu = woGu CONCAT $woUQ
            set woGl = woGl CONCAT $woLQ
            set randGm = randGm CONCAT $randMean
            set randGs = randGs CONCAT $randSigma
            set randGu = randGu CONCAT $randUQ
            set randGl = randGl CONCAT $randLQ
           


PNeErrorg
    expand 1.0001
    define TeX_strings 1
    #graphin' time (tm)
    #this is run once for each image (3,3m,4,4m,etc)
    #can now use vectors defined above in PNeErrorc to graph things
   
    #data in vectors is one row / range value, ranges are 10-60, by 5s   
    #will graph in terms of arcseconds, not pixels
    set rangesG = ra*1.45
   
   
    device postencap '$!v.0.med$imag.we.ps'
#    device x11
    erase
    #top window left: 1+ PNe Medians
        window 2 -3 1 3
        ctype black
        vecminmax wGm wmin wmax
        lim 10 90 $($wmin-$($wmin*.5)) $($wmax+$($wmax*.5))
        box
        connect rangesG wGm
#        #temp for initial testing
#        relocate $(rangesG[0]) $(wGm[0])
#        #dot $($wmin-$($wmin*.5)) $($wmax+$($wmax*.5))
        ylabel Random w/ 1+ PNe
        #sigma error bars
        ctype blue
        set sig = abs(wGs)
        errorbar rangesG wGm sig 2
        errorbar rangesG wGm sig 4
        #U/L quartile error bars
        ctype red
        set R2u = abs(wGu - wGm)
        set R2l = abs(wGl - wGm)
        errorbar rangesG wGm R2u 2
        errorbar rangesG wGm R2l 4
    #basic labels
    lim 10 90 0 10.2
    ctype black
    relocate 25 11.6
    label $imag  - 
    ctype red
    relocate 110 11.6
    label  red-upper/lower quartile
    ctype blue
    relocate 70 11.6
    label blue-sigma
    ctype black
    relocate 25 10.5
    label average medians
    relocate 143 10.5
    label delta plots
    relocate 87 10.3
    label \it{$(dimen(randMt))+ iterations}
       
    #top middle window: Random PNe Medians
        window 2 -3 1 2
        ctype black
        vecminmax randGm randmin randmax
        lim 10 90 $($randmin-$($randmin*.5)) $($randmax+$($randmax*.5))
        box
        connect rangesG randGm
#        #temp for initial testing
#        relocate $(rangesG[0]) $(randGm[0])
#        dot
        ylabel Random
        #sigma error bars
        ctype blue
        set sig = abs(randGs)
        errorbar rangesG randGm sig 2
        errorbar rangesG randGm sig 4
        #U/L quartile error bars
        ctype red
        set Ru = abs(randGu - randGm)
        set Rl = abs(randGl - randGm)
        errorbar rangesG randGm Ru 2
        errorbar rangesG randGm Rl 4
   
    #bottom middle window: 0 PNe Medians
        window 2 -3 1 1
        ctype black
        vecminmax woGm womin womax
        lim 10 90 $($womin-$($womin*.5)) $($womax+$($womax*.5))
        box
        connect rangesG woGm
#        #temp for initial testing
#        relocate $(rangesG[0]) $(woGm[0])
#        dot
        ylabel Random w/ 0 PNe
        xlabel half-box width (arcseconds)
        #sigma error bars
        ctype blue
        set sig = abs(woGs)
        errorbar rangesG woGm sig 2
        errorbar rangesG woGm sig 4
        #U/L quartile error bars
        ctype red
        set wou = abs(woGu-woGm)
        set wol = abs(woGl-woGm)
        errorbar rangesG woGm wou 2
        errorbar rangesG woGm wol 4

   
    #do the deltas!
    #calculate deltas
    set Rw = wGm - randGm
    set Rwo = woGm - randGm
    set Rwwo = wGm - woGm
    #find limits
    vecminmax Rw wmin wmax
    vecminmax Rwo womin womax
    vecminmax Rwwo wwomin wwomax
    set dimen(chk) = 6
    set chk[0]=$wmin
    set chk[1]=$wmax
    set chk[2]=$womin
    set chk[3]=$womax
    set chk[4]=$wwomin
    set chk[5]=$wwomax
    vecminmax chk ymin ymax
    #if ($twomin < $onmin && $twomin < $offmin) { define ymin $twomin }
    #if ($onmin < $twomin && $onmin < $offmin) { define ymin $onmin }
    #if ($offmin < $twomin && $offmin < $onmin) { define ymin $offmin }
    #if ($twomax > $onmax && $twomax > $offmax) { define ymax $twomax }
    #if ($onmax > $twomax && $onmax > $offmax) { define ymax $onmax }
    #if ($offmax > $twomax && $offmax > $onmax) { define ymax $offmax }
    #add some buffer around the edges for error bars
    define ymax $($ymax + $(abs($ymax-$ymin)/2))
    define ymin $($ymin - $(abs($ymax-$ymin)/2))
   
    #deltas on the right side
    #top window : delta (R - 1)
        window 2 -3 2 3
        ctype black
        lim 10 90 $ymin $ymax
        box
        connect rangesG Rw
        ylabel Delta(1 - R)
        #sigma error bars
        ctype blue
        set sig = abs(abs(randGs) + abs(wGs))
        errorbar rangesG Rw sig 2
        errorbar rangesG Rw sig 4
        #U/L quartile error bars
        ctype red
        set u = abs(abs(randGu-randGm)+abs(wGu-wGm))
        set l = abs(abs(randGl-randGm)+abs(wGl-wGm))
        errorbar rangesG Rw u 2
        errorbar rangesG Rw l 4
   
    #middle window : delta (R - 0)
        window 2 -3 2 2
        ctype black
        lim 10 90 $ymin $ymax
        box
        connect rangesG Rwo
        ylabel Delta(0 - R)
        #sigma error bars
        ctype blue
        set sig = abs(abs(randGs)+abs(woGs))
        errorbar rangesG Rwo sig 2
        errorbar rangesG Rwo sig 4
        #U/L quartile error bars
        ctype red
        set u = abs(abs(randGu - randGm)+abs(woGu-woGm))
        set l = abs(abs(randGl - randGm)+abs(woGl-woGm))
        errorbar rangesG Rwo u 2
        errorbar rangesG Rwo l 4
       
    #bottom window : delta (1 - 0)
        window 2 -3 2 1
        ctype black
        lim 10 90 $ymin $ymax
        box
        connect rangesG Rwwo
        ylabel Delta(1 - 0)
        xlabel half-box width (arcseconds)
        #sigma error bars
        ctype blue
        set sig = abs(abs(wGs)+abs(woGs))
        errorbar rangesG Rwwo sig 2
        errorbar rangesG Rwwo sig 4
        #U/L quartile error bars
        ctype red
        set u = abs(abs(wGu - wGm)+abs(woGu-woGm))
        set l = abs(abs(wGl - wGm)+abs(woGl-woGm))
        errorbar rangesG Rwwo u 2
        errorbar rangesG Rwwo l 4
       

qc
    !convert 8.0.med3.fits.we.ps 8.0.med3.fits.we.ps.gif
    !convert 8.0.med3m.fits.we.ps 8.0.med3m.fits.we.ps.gif
    !convert 8.0.med4.fits.we.ps 8.0.med4.fits.we.ps.gif
    !convert 8.0.med4m.fits.we.ps 8.0.med4m.fits.we.ps.gif
    !convert 8.0.med7.fits.we.ps 8.0.med7.fits.we.ps.gif
    !convert 8.0.med7m.fits.we.ps 8.0.med7m.fits.we.ps.gif
    !convert 8.0.medSub.fits.we.ps 8.0.medSub.fits.we.ps.gif
    !convert 8.0.medSubm.fits.we.ps 8.0.medSubm.fits.we.ps.gif
    !cp *.gif /astroweb_data/web/steven/PNe8/
    !cp *.gif ~/Virgo/Web/PNe8/