/packmule_raid/steven/PNe15/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])
define itt 21
set it = 1,$itt
#set ra = 10,60,5
set ra = { 1 3 5 7 9 11 13 15 25 35 45 55 65 }
set cs = 0,10,2
#temp for testing
# set im = { 3.fits }
# set ra = { 1 }
# set cs = { 2 4 6 8 10 }
define ymaxset 25
define yminset -10
#separate by contamination level
foreach c cs {
#separate main runings, next by
imagee
foreach imag im {
if ('$imag' ==
'"3.fits"' || '$imag' == '"3m.fits"') {
define ymaxset 40
define yminset -15
}
if ('$imag' ==
'"4.fits"' || '$imag' == '"4m.fits"') {
define ymaxset 4
define yminset -6
}
if ('$imag' ==
'"7.fits"') {
define ymaxset 1000
define yminset -50
}
if ('$imag' ==
'"7m.fits"') {
define ymaxset 40
define yminset -50
}
if ('$imag' ==
'"Sub.fits"') {
define ymaxset 1000
define yminset -100
}
if ('$imag' ==
'"Subm.fits"') {
define ymaxset 50
define yminset -30
}
#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
#add on for
delta support
set
dimen(DwGm) = 0
set
dimen(DwGs) = 0
set
dimen(DwGu) = 0
set
dimen(DwGl) = 0
set
dimen(DwoGm) = 0
set
dimen(DwoGs) = 0
set
dimen(DwoGu) = 0
set
dimen(DwoGl) = 0
set
dimen(DwwoGm) = 0
set
dimen(DwwoGs) = 0
set
dimen(DwwoGu) = 0
set
dimen(DwwoGl) = 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
}
#end of contaminations 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.C3.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
# C3 is contam, 0-10, by 1s
# file will look like:
# write
"$!v.$!itn.Im$!imag.R$!rang.C$!c.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.C$!c.med)) {
data $v.$itn.Im$imag.R$rang.C$!c.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
#compute
deltas for this range
set Dw = wMt -
randMt
set Dwo = woMt
- randMt
set Dwwo = wMt
- woMt
#finding sigma:
stats wMt
wMean wSigma wK
stats woMt
woMean woSigma woK
stats randMt
randMean randSigma randK
stats Dw
DwMean DwSigma DwK
stats Dwo
DwoMean DwoSigma DwoK
stats Dwwo
DwwoMean DwwoSigma DwwoK
#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)])
sort {Dw}
define DwMed
(Dw[DIMEN(Dw)/2])
define DwLQ
(Dw[0.25*DIMEN(Dw)])
define DwUQ
(Dw[0.75*DIMEN(Dw)])
sort {Dwo}
define DwoMed
(Dwo[DIMEN(Dwo)/2])
define DwoLQ
(Dwo[0.25*DIMEN(Dwo)])
define DwoUQ
(Dwo[0.75*DIMEN(Dwo)])
sort {Dwwo}
define DwwoMed
(Dwwo[DIMEN(Dwwo)/2])
define DwwoLQ
(Dwwo[0.25*DIMEN(Dwwo)])
define DwwoUQ
(Dwwo[0.75*DIMEN(Dwwo)])
#write results
to data file
!rm
"$!v.0.Im$!imag.R$!rang.C$!c.stats"
#writerand1Gm
header
write
"$!v.0.Im$!imag.R$!rang.C$!c.stats" "#wPNe woPNe rand Dw Dwo Dwoo -
data rows go: (1)Mean, (2)Sigma, (3)Kurtosis, (4)Median, (5)Lower
Quartile, (6)Upper Quartile"
write
"$!v.0.Im$!imag.R$!rang.C$!c.stats" $wMean $woMean $randMean $DwMean
$DwoMean $DwwoMean
write
"$!v.0.Im$!imag.R$!rang.C$!c.stats" $wSigma $woSigma $randSigma
$DwSigma $DwoSigma $DwwoSigma
write
"$!v.0.Im$!imag.R$!rang.C$!c.stats" $wK $woK $randK $DwK $DwoK $DwwoK
write
"$!v.0.Im$!imag.R$!rang.C$!c.stats" $wMed $woMed $randMed $DwMed
$DwoMed $DwwoMed
write
"$!v.0.Im$!imag.R$!rang.C$!c.stats" $wLQ $woLQ $randLQ $DwLQ $DwoLQ
$DwwoLQ
write
"$!v.0.Im$!imag.R$!rang.C$!c.stats" $wUQ $woUQ $randUQ $DwUQ $DwoUQ
$DwwoUQ
#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
#add'l delta
support
set DwGm =
DwGm CONCAT $DwMean
set DwGs =
DwGs CONCAT $DwSigma
set DwGu =
DwGu CONCAT $DwUQ
set DwGl =
DwGl CONCAT $DwLQ
set DwoGm =
DwoGm CONCAT $DwoMean
set DwoGs =
DwoGs CONCAT $DwoSigma
set DwoGu =
DwoGu CONCAT $DwoUQ
set DwoGl =
DwoGl CONCAT $DwoLQ
set DwwoGm =
DwwoGm CONCAT $DwwoMean
set DwwoGs =
DwwoGs CONCAT $DwwoSigma
set DwwoGu =
DwwoGu CONCAT $DwwoUQ
set DwwoGl =
DwwoGl CONCAT $DwwoLQ
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.C$!c.we.ps'
# device x11
erase
#top window left: 1+ PNe Medians
window 2 -3 1 3
ctype black
vecminmax wGm wmin wmax
lim -1 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 91 10.3
label \it{$(dimen(randMt))+ trials}
expand 1.5
relocate 97 9.1
label $($(int($c))*10)%
expand 1.0001
relocate 92 8.4
label contam.
#top middle window: Random PNe Medians
window 2 -3 1 2
ctype black
vecminmax randGm randmin randmax
lim -1 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 -1 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 - ANTIQUATED
##set Rw = wGm - randGm
##set Rwo = woGm - randGm
##set Rwwo = wGm - woGm
#find limits
vecminmax DwGm wmin wmax
vecminmax DwoGm womin womax
vecminmax DwwoGm 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))
define ymax $ymaxset
define ymin $yminset
#deltas on the right side
#top window : delta (R - 1)
window 2 -3 2 3
ctype black
lim -1 90 $ymin $ymax
box
connect rangesG DwGm
ylabel Delta(1 - R)
#sigma error bars
ctype blue
set sig = abs(DwGs)
errorbar rangesG DwGm sig 2
errorbar rangesG DwGm sig 4
#U/L quartile error bars
ctype red
set u = abs(DwGu-DwGm)
set l = abs(DwGl-DwGm)
errorbar rangesG DwGm u 2
errorbar rangesG DwGm l 4
#middle window : delta (R - 0)
window 2 -3 2 2
ctype black
lim -1 90 $ymin $ymax
box
connect rangesG DwoGm
ylabel Delta(0 - R)
#sigma error bars
ctype blue
set sig = abs(DwoGs)
errorbar rangesG DwoGm sig 2
errorbar rangesG DwoGm sig 4
#U/L quartile error bars
ctype red
set u = abs(DwoGu - DwoGm)
set l = abs(DwoGl - DwoGm)
errorbar rangesG DwoGm u 2
errorbar rangesG DwoGm l 4
#bottom window : delta (1 - 0)
window 2 -3 2 1
ctype black
lim -1 90 $ymin $ymax
box
connect rangesG DwwoGm
ylabel Delta(1 - 0)
xlabel half-box width (arcseconds)
#sigma error bars
ctype blue
set sig = abs(DwwoGs)
errorbar rangesG DwwoGm sig 2
errorbar rangesG DwwoGm sig 4
#U/L quartile error bars
ctype red
set u = abs(DwwoGu - DwwoGm)
set l = abs(DwwoGl - DwwoGm)
errorbar rangesG DwwoGm u 2
errorbar rangesG DwwoGm l 4
qc
echo converting...
foreach c cs {
!convert
15.0.med3.fits.C$!c.we.ps 15.0.med3.fits.C$!c.we.ps.gif
!convert
15.0.med3m.fits.C$!c.we.ps 15.0.med3m.fits.C$!c.we.ps.gif
!convert
15.0.med4.fits.C$!c.we.ps 15.0.med4.fits.C$!c.we.ps.gif
!convert
15.0.med4m.fits.C$!c.we.ps 15.0.med4m.fits.C$!c.we.ps.gif
!convert
15.0.med7.fits.C$!c.we.ps 15.0.med7.fits.C$!c.we.ps.gif
!convert
15.0.med7m.fits.C$!c.we.ps 15.0.med7m.fits.C$!c.we.ps.gif
!convert
15.0.medSub.fits.C$!c.we.ps 15.0.medSub.fits.C$!c.we.ps.gif
!convert
15.0.medSubm.fits.C$!c.we.ps 15.0.medSubm.fits.C$!c.we.ps.gif
}
echo ...copying...
!cp *.gif /astroweb_data/web/steven/PNe15/
!cp *.gif ~/Virgo/Web/PNe15/
echo ...done