# This is a beautiful day and it will stay just as lovely setup #include "/home/znadia1/sm/addons.sm" play include "/data1/lib/sm/spectro.sm" play #include "/home/znadia1/sm/multiple_matching.sm" play eightwin 1 # eightwin windownumber # set up for eight rectangular windows on laser_p if ($1==1) {location 3000 15000 25000 32050} if ($1==2) {location 3000 15000 17250 24300} if ($1==3) {location 3000 15000 9500 16550} if ($1==4) {location 3000 15000 1750 8800} if ($1==5) {location 19500 31500 25000 32050} if ($1==6) {location 19500 31500 17250 24300} if ($1==7) {location 19500 31500 9500 16550} if ($1==8) {location 19500 31500 1750 8800} eightwin2 1 # eightwin windownumber # like eightwin but tighter, no space for right side labels if ($1==1) {location 3000 14500 25000 32050} if ($1==2) {location 3000 14500 17250 24300} if ($1==3) {location 3000 14500 9500 16550} if ($1==4) {location 3000 14500 1750 8800} if ($1==5) {location 15500 27000 25000 32050} if ($1==6) {location 15500 27000 17250 24300} if ($1==7) {location 15500 27000 9500 16550} if ($1==8) {location 15500 27000 1750 8800} eightwin3 1 # eightwin windownumber # set up for eight rectangular windows on laser_p if ($1==1) {location 3500 15500 24200 31050} if ($1==2) {location 3500 15500 16650 23500} if ($1==3) {location 3500 15500 9100 15950} if ($1==4) {location 3500 15500 1550 8400} if ($1==5) {location 19500 31500 24200 31050} if ($1==6) {location 19500 31500 16650 23500} if ($1==7) {location 19500 31500 9100 15950} if ($1==8) {location 19500 31500 1550 8400} eightwin4 1 # eightwin windownumber # set up for eight rectangular windows on laser_p # enough room for x labels if ($1==1) {location 3500 15500 25058 31050} if ($1==2) {location 3500 15500 17222 23214} if ($1==3) {location 3500 15500 9386 15378} if ($1==4) {location 3500 15500 1550 7542} if ($1==5) {location 19500 31500 25058 31050} if ($1==6) {location 19500 31500 17222 23214} if ($1==7) {location 19500 31500 9386 15378} if ($1==8) {location 19500 31500 1550 7542} eightlabel3 2 # eightlabel3 y label # header label for eightwin3 set up, centered at y/10 location 3500 31500 24200 31050 limits 0 10 0 10 relocate 5 $1 putlabel 5 $2 figlabel 1 # figlabel n # put label "Figure n" in upper right corner of laser_p plot location 0 32768 0 32768 limits 0 32768 0 32768 relocate 31500 32000 expand 0.95 putlabel 4 {\ti Fig. $1} figlabel1 1 # figlabel1 n # Figure label for "twowin 1" set-up location 0 32768 0 32768 limits 0 32768 0 32768 relocate 28500 31500 expand 0.95 putlabel 4 {\ti Fig. $1} figlabel2 1 # figlabel2 n # Figure label for "twowin 2" set-up location 0 32768 0 32768 limits 0 32768 0 32768 relocate 28500 17000 expand 0.95 putlabel 4 {\ti Fig. $1} fileavg 8 # fileavg name1 name2 nfiles xcol ycol x y sig # average results from multiple files set y=0 set sig=0 da $11$2 read <$6 $4 $7 $5> set $8 = $7**2 do i=2,$3 { da $1"$!i"$2 read _y $5 set $7 = $7 + _y set $8 = $8 + _y**2 } set $7 = $7/$3 set $8 = sqrt($8/$3-$7**2) fourwin 1 # fourwin windownumber # set up for four large square windows on laser_p # like sixwin2, but slightly more space between levels, # allowing axis labels for upper panels if ($1==1) {location 5000 16380 21500 30000} if ($1==2) {location 18950 30330 21500 30000} if ($1==3) {location 5000 16380 10500 19000} if ($1==4) {location 18950 30330 10500 19000} head 3 # set one vector equal to the head of another # head n vec1 vec2 # on output, vec1 contains the n first elements of vec2 set dimen($2)=$1 do _i=0,$1-1 {$2[_$i]=$3[_$i]} laser # laser printer, square mode device postscript laser_p # laser printer, portrait mode device postport laser_l # laser printer, landscape mode device postland laser_s # laser printer, square mode device postscript laser_g # laser printer, for things in grey (square) device postgray laser_lp 1 # postscript file, laser landscape mode # laser_lp filename device postlandfile $1 laser_pp 1 # postscript file, laser portrait mode # laser_pp filename device postportfile $1 laser_sp 1 # postscript file, laser square mode # laser_sp filename device postencap $1 laser_gp 1 # postscript file, for things in grey (square) device postgrayfile $1 lbug # temporary fix for lweight bug overload lweight 1 \n macro lweight 1 {define lweight $1 \n LWEIGHT $1} ninelabel 2 # ninelabel y label # header label for ninewin set up, centered at y/10 location 3000 31000 23280 30000 limits 0 10 0 10 relocate 5 $1 putlabel 5 $2 ninewin 1 # ninewin windownumber # set up for nine square windows on laser_p if ($1==1) {location 3000 12000 23280 30000} if ($1==2) {location 3000 12000 15280 22000} if ($1==3) {location 3000 12000 7280 14000} if ($1==4) {location 12500 21500 23280 30000} if ($1==5) {location 12500 21500 15280 22000} if ($1==6) {location 12500 21500 7280 14000} if ($1==7) {location 22000 31000 23280 30000} if ($1==8) {location 22000 31000 15280 22000} if ($1==9) {location 22000 31000 7280 14000} orlabel 6 # orlabel y pexpand nsides istyle lexpand label # put label and point outside upper right corner, at y=y foreach var1 {gx1 gx2 gy1 gy2} {define $var1 |} define ngx2 ($gx1 + 1.2*($gx2-$gx1)) location $gx1 $ngx2 $gy1 $gy2 relocate 9 $1 expand $2 ptype $3 $4 dot relocate 9.3 $1 expand $5 putlabel 6 $6 location $gx1 $gx2 $gy1 $gy2 foreach var1 {gx1 gx2 gy1 gy2 ngx2} {define $var1 delete} define var1 delete oralabel 7 # oralabel y pexpand pangle nsides istyle lexpand label # put label and angled point outside upper right corner, at y=y foreach var1 {gx1 gx2 gy1 gy2} {define $var1 |} define ngx2 ($gx1 + 1.2*($gx2-$gx1)) location $gx1 $ngx2 $gy1 $gy2 relocate 9 $1 expand $2 angle $3 ptype $4 $5 dot angle 0 relocate 9.3 $1 expand $6 putlabel 6 $7 location $gx1 $gx2 $gy1 $gy2 foreach var1 {gx1 gx2 gy1 gy2 ngx2} {define $var1 delete} define var1 delete per_s # pericom call, with location for square plot pericom location 3500 24447 3500 31000 putdate # put date in lower right hand corner location 0 32768 0 32768 limits 0 32768 0 32768 expand 0.55 relocate 32600 400 putlabel 4 {\ti $date} radtick 6 # put tick marks on a radial axis # radtick dtick rmax ltick angle xdir ydir define _cos (cos($4*pi/180.)) define _sin (sin($4*pi/180.)) do _rad=$1,$2,$1 { rel ($_rad*$_cos) ($_rad*$_sin) dra ($_rad*$_cos+$5*$3*$_sin) ($_rad*$_sin+$6*$3*$_cos) } define _cos delete define _sin delete define _rad delete ran_num 14 # generate quick-and-dirty uniform deviates # random seed nran x [newseed] # returns a vector of nran random numbers in x[] # employs linear congruential pseudo-random generator with # starting seed given by first argument # if fourth argument exists, will return value for new seed, # which can be employed in a later call define _im 7875 define _ia 421 define _ic 1663 set dimen($3)=$2 define _jran ($1*$_ia+$_ic - $_im*int(($1*$_ia+$_ic)/$_im)) do _iran=0,$2-1 { set $3[$_iran]=$_jran/$_im define _jran ($_jran*$_ia+$_ic - $_im*int(($_jran*$_ia+$_ic)/$_im)) } if ($?4) {define $4 $_jran} foreach v (_im _ia _ic _jran _iran) {define $v delete} sixlabel 2 # sixlabel y label # header label for sixwin setup, centered at y/10 location 5000 29000 23100 30000 limits 0 10 0 10 relocate 5 $1 putlabel 5 $2 sixwin 1 # sixwin windownumber # set up for six square windows on laser_p if ($1==1) {location 5000 15000 22530 30000} if ($1==2) {location 19000 29000 22530 30000} if ($1==3) {location 5000 15000 12860 20330} if ($1==4) {location 19000 29000 12860 20330} if ($1==5) {location 5000 15000 3190 10660} if ($1==6) {location 19000 29000 3190 10660} sixwin2 1 # sixwin2 windownumber # set up for six large square windows on laser_p if ($1==1) {location 5000 16380 21500 30000} if ($1==2) {location 18950 30330 21500 30000} if ($1==3) {location 5000 16380 11500 20000} if ($1==4) {location 18950 30330 11500 20000} if ($1==5) {location 5000 16380 1500 10000} if ($1==6) {location 18950 30330 1500 10000} sixlabel2 2 # sixlabel y label # header label for sixwin2 setup, centered at y/10 location 5000 30330 21500 30000 limits 0 10 0 10 relocate 5 $1 putlabel 5 $2 sixwin3 1 # sixwin3 windownumber # set up for six large square windows on laser_p, with space # for axis labels below bottom panels if ($1==1) {location 5000 15710 22000 30000} if ($1==2) {location 18300 29010 22000 30000} if ($1==3) {location 5000 15710 12500 20500} if ($1==4) {location 18300 29010 12500 20500} if ($1==5) {location 5000 15710 3000 11000} if ($1==6) {location 18300 29010 3000 11000} sixlabel3 2 # sixlabel y label # header label for sixwin2 setup, centered at y/10 location 5000 29010 22000 30000 limits 0 10 0 10 relocate 5 $1 putlabel 5 $2 sixwin4 1 # sixwin4 windownumber # set up for six rectangular windows on laser_p if ($1==1) {location 4000 16000 23000 30050} if ($1==2) {location 19500 31500 23000 30050} if ($1==3) {location 4000 16000 14250 21300} if ($1==4) {location 19500 31500 14250 21300} if ($1==5) {location 4000 16000 5500 12550} if ($1==6) {location 19500 31500 5500 12550} sixlabel4 2 # sixlabel4 y label # header label for sixwin4 setup, centered at y/10 location 4000 31500 23000 30050 limits 0 10 0 10 relocate 5 $1 putlabel 5 $2 sixwinl 1 # sixwinl windownumber # set up for six windows on laser_l, rotation of sixwin2 # first row is 1,2,3 ; second row is 4,5,6 if ($1==1) {location 2500 11000 16950 28330} if ($1==2) {location 12500 21000 16950 28330} if ($1==3) {location 22500 31000 16950 28330} if ($1==4) {location 2500 11000 3000 14380} if ($1==5) {location 12500 21000 3000 14380} if ($1==6) {location 22500 31000 3000 14380} spairs 7 # connect pairs in a slice # pairs x1 y1 x2 y2 z zmin zmax. connect (x1,y1) to (x2,y2) if # zmin < z < zmax DEFINE 9 { fx1 fx2 fy1 fy2 gx1 gx2 gy1 gy2 ptype angle expand aspect } FOREACH 8 { $!!9 } { DEFINE $8 | } ASPECT 1 SET _dx$0=($3 - $1)*($gx2 - $gx1)/($fx2 - $fx1) SET _dy$0=($4 - $2)*($gy2 - $gy1)/($fy2 - $fy1) PTYPE 2 0 SET _a$0=(_dx$0 == 0 ? (_dy$0 > 0 ? PI/2 : -PI/2) : \ (_dy$0 > 0 ? ATAN(_dy$0/_dx$0) : ATAN(_dy$0/_dx$0) + PI)) ANGLE 180/pi*_a$0 EXPAND SQRT(_dx$0**2 + _dy$0**2)/(2*128) POINTS (($1 + $3)/2) (($2 + $4)/2) if ($5>$6 && $5<$7) FOREACH 8 { angle aspect expand ptype } { $8 $$8 } FOREACH 8 { $!!9 } { DEFINE $8 DELETE } FOREACH 8 ( _a _dx _dy ) { DELETE $8$0 } squarelim 2 # squarelim x y # like "limits x y" but with equal axis lengths limits $1 $2 define _xr ($fx2-$fx1) define _yr ($fy2-$fy1) if ($_xr>$_yr) { define _min (($fy2+$fy1-$_xr)/2.) define _max (($fy2+$fy1+$_xr)/2.) limits $1 $_min $_max } else { define _min (($fx2+$fx1-$_yr)/2.) define _max (($fx2+$fx1+$_yr)/2.) limits $_min $_max $2 } define _xr delete define _yr delete define _min delete define _max delete sun ## sunview window, in convenient location device sunview -Ws 650 525 -Wp 450 325 sun_l ## sunview window, same dimensions as landscape postscript device sunview -Ws 837 625 -Wp 300 250 sun_p ## sunview window, same dimensions as portrait postscript device sunview -Ws 625 837 -Wp 500 50 sun_s ## square sun window device sunview -Ws 525 525 -Wp 575 325 endpiece 3 # set one vector equal to the tail of another # endpiece n vec1 vec2 # on output, vec1 contains the n last elements of vec2 set dimen($2)=$1 define _dim (dimen($3)) do _i=0,$1-1 {set $2[$_i]=$3[($_dim - $1 + $_i)]} define _i delete define _dim delete thin 3 # create a "thinned" version of a vector for plotting points # thin x n y # on return, vector y contains every n'th element of x define _n1 (int((dimen($1)-1)/$2)+1) set dimen($3) = $_n1 do _i=0,$_n1-1 {set $3[$_i]=$1[($2*$_i)]} define _n1 delete define _i delete thin2 4 # "thin" the latter portion of a vector for plotting points # thin2 x m n y # on return, vector y contains the first m elements of x and # every n'th element of x thereafter define _n1 ($2+int((dimen($1)-$2)/$3)) set dimen($4) = $_n1 do _i=0,$2-1 {set $4[$_i]=$1[$_i]} do _i=$2,$_n1-1 {set $4[$_i]=$1[($2-1+$3*($_i-$2+1))]} define _n1 delete define _i delete thin3 4 # "thin" the central portion of a vector for plotting points # thin3 x m n y # on return, vector y contains the first m elements of x, the # last m elements of x, and every n'th element of x in between define _n1 (2*$2+int((dimen($1)-2*$2)/$3)) set dimen($4) = $_n1 do _i=0,$2-1 {set $4[$_i]=$1[$_i]} do _i=$2,$_n1-1-$2 {set $4[$_i]=$1[($2-1+$3*($_i-$2+1))]} do _i=$_n1-$2,$_n1-1 {set $4[$_i]=$1[(dimen($1)+$_i-$_n1)]} define _n1 delete define _i delete twowin 1 # twowin windownumber # set up for two vertically stacked windows on laser_p if ($1==1) {location 8000 25000 19500 30000} if ($1==2) {location 8000 25000 5000 15500} twowin1 1 # twowin1 windownumber # More stretched out if ($1==1) {location 3000 30000 17500 30000} if ($1==2) {location 3000 30000 3000 15500} ullabel 6 # ullabel y pexpand nsides istyle lexpand label # put label and point in upper left corner, at y=y relocate 1 $1 expand $2 ptype $3 $4 dot relocate 1.3 $1 expand $5 putlabel 6 $6 urlabel 6 # urlabel y pexpand nsides istyle lexpand label # put label and point in upper right corner, at y=y relocate 9 $1 expand $2 ptype $3 $4 dot relocate 8.7 $1 expand $5 putlabel 4 $6 uralabel 7 # uralabel y pexpand pangle nsides istyle lexpand label # put label and angled point in upper right corner, at y=y relocate 9 $1 expand $2 angle $3 ptype $4 $5 dot angle 0 relocate 8.7 $1 expand $6 putlabel 4 $7 ullinlbl 3 # ullinlbl y ltype label # put label and line in upper left corner, at y=y relocate 0.5 $1 ltype $2 draw 1.5 $1 ltype 0 relocate 1.8 $1 putlabel 6 $3 urlinlbl 3 # urlinlbl y ltype label # put label and line in upper right corner, at y=y relocate 9.5 $1 ltype $2 draw 8.5 $1 ltype 0 relocate 8.2 $1 putlabel 4 $3 vfield2 4 # plot a vector field: vfield x y len angle # len is length of vector (in x-axis user units) # angle is in +ve direction from +x axis # arrows start at (x,y) DEFINE ptype | PTYPE { 600 0 m 450 100 600 0 450 -100 } DEFINE angle | ANGLE $4 foreach v {gx1 gx2 fx1 fx2} {DEFINE $v |} DEFINE _cfac (($gx2-$gx1)/(600.*($fx2-$fx1))) DEFINE expand | EXPAND ($_cfac*$3) POINTS $1 $2 FOREACH v { angle expand ptype } {$v $$v DEFINE $v DELETE} DEFINE _cfac DELETE vfield3 4 # plot a vector field: vfield x y len angle # len is length of vector (in y-axis user units) # angle is in +ve direction from +x axis # arrows start at (x,y) DEFINE ptype | PTYPE { 600 0 m 450 100 600 0 450 -100 } DEFINE angle | ANGLE $4 foreach v {gy1 gy2 fy1 fy2} {DEFINE $v |} DEFINE _cfac (($gy2-$gy1)/(600.*($fy2-$fy1))) DEFINE expand | EXPAND ($_cfac*$3) POINTS $1 $2 FOREACH v { angle expand ptype } {$v $$v DEFINE $v DELETE} DEFINE _cfac DELETE xterm ## xterm window, in convenient location device x11 -bg white -geometry 650x525+450+325 xterm_b #A really big one device x11 -bg white -geometry 1200x800 xterm_l ## xterm window, same dimensions as landscape postscript device x11 -bg white -geometry 837x625+300+230 xterm_p ## xterm window, same dimensions as portrait postscript device x11 -bg white -geometry 625x837+500+30 xterm_h ## xterm window, square ## shrunk down for home device x11 -bg white -geometry 500x500+450+30 xterm_s ## square xterm window device x11 -bg white -geometry 800x800+325+50 xterm_ss ## small square xterm window device x11 -bg white -geometry 525x525+525+50 display 1 #Switch between display windows device x11 -dev $1 12win 1 # 12win windownumber # set up for 12 square windows (4x3) on laser_l # window 1 is upper left, 2 is immediately below, etc. if ($1==1) {location 2000 8720 21500 30500} if ($1==2) {location 2000 8720 11500 20500} if ($1==3) {location 2000 8720 1500 10500} if ($1==4) {location 9500 16220 21500 30500} if ($1==5) {location 9500 16220 11500 20500} if ($1==6) {location 9500 16220 1500 10500} if ($1==7) {location 17000 23720 21500 30500} if ($1==8) {location 17000 23720 11500 20500} if ($1==9) {location 17000 23720 1500 10500} if ($1==10) {location 24500 31220 21500 30500} if ($1==11) {location 24500 31220 11500 20500} if ($1==12) {location 24500 31220 1500 10500} 16win 1 # 16win windownumber # set up for 16 square windows (4x4) on laser_p # window 1 is upper left, 2 is immediately below, etc. if ($1==1) {location 2000 9100 24770 30000 } if ($1==2) {location 2000 9100 19020 24250 } if ($1==3) {location 2000 9100 13270 18500 } if ($1==4) {location 2000 9100 7520 12750 } if ($1==5) {location 9750 16850 24770 30000} if ($1==6) {location 9750 16850 19020 24250} if ($1==7) {location 9750 16850 13270 18500} if ($1==8) {location 9750 16850 7520 12750} if ($1==9) {location 17500 24600 24770 30000} if ($1==10) {location 17500 24600 19020 24250} if ($1==11) {location 17500 24600 13270 18500} if ($1==12) {location 17500 24600 7520 12750} if ($1==13) {location 25250 32350 24770 30000} if ($1==14) {location 25250 32350 19020 24250} if ($1==15) {location 25250 32350 13270 18500} if ($1==16) {location 25250 32350 7520 12750} 16label 2 # sixlabel y label # header label for 16win setup, centered at y/10 location 2000 32350 24770 30000 limits 0 10 0 10 relocate 5 $1 putlabel 5 $2 histshed 5 # shade region between x y and x2 y2 with n lines # syntax: shed x y x2 y2 n SET _x = reverse($3) SET _x = $1 CONCAT _x SET _y = reverse($4) SET _y = $2 CONCAT _y DEFINE _n (INT(32767/$5)) SHADE HISTOGRAM $_n _x _y FOREACH v ( _x _y ) { DELETE $v } DEFINE _n DELETE readcoor 2 #Read in ra and dec read '%2d%2d%2d.%d%1s%2d%2d%2d' \ {_hour _ramin _rasec _fracsec _sign _deg _decmin _decsec} set $1 = _hour + _ramin/60. + _rasec/3600. + _fracsec/36000. set $2 = _deg + _decmin/60. + _decsec/3600. set $2 = ((_sign == '-')? (-$2) : ($2)) delete _hour delete _ramin delete _rasec delete _fracsec delete _sign delete _deg delete _decmin delete _decsec oreadcoor 2 #Read in ra and dec read '%*10[]%2d%2d%2d.%d%1s%2d%2d%2d' \ {_hour _ramin _rasec _fracsec _sign _deg _decmin _decsec} set $1 = _hour + _ramin/60. + _rasec/3600. + _fracsec/36000. set $2 = _deg + _decmin/60. + _decsec/3600. set $2 = ((_sign == '-')? (-$2) : ($2)) delete _hour delete _ramin delete _rasec delete _fracsec delete _sign delete _deg delete _decmin delete _decsec zreadcoor 4 #Read in ra, dec, mag, and redshift. For ZCAT format read '%*11[]%8s%7s%5s%5s' {_rastring _decstring _magstring _czstring} set _hour = substr(_rastring,0,2) set _ramin = substr(_rastring,2,2) set _rasec = substr(_rastring,4,4) set _sign = substr(_decstring,0,1) set _deg = substr(_decstring,1,2) set _decmin = substr(_decstring,3,2) set _decsec = substr(_decstring,5,2) set _mag = (_magstring == ' ')? ('-999') : (_magstring) set _cz = (_czstring == ' ')? ('-999') : (_czstring) set _decsec = (_decsec == ' ')?('00'):(_decsec) foreach i {_hour _ramin _rasec _deg _decmin _decsec} { set _$i = atof($i)} set $1 = __hour + __ramin/60. + __rasec/3600. set $2 = __deg + __decmin/60. + __decsec/3600. set $2 = ((_sign == '-')? (-$2) : ($2)) set $3 = atof(_mag) set $4 = atof(_cz) foreach i {_hour _ramin _rasec _deg _decmin _decsec} { delete $i delete _$i} delete _rastring delete _decstring delete _magstring delete _czstring delete _mag delete _cz readvel 2 #Read in redshift and flux read '%*24[]%4d%*32[]%5d' {_temp1 _temp2} set $1 = 10.**(_temp1/1000. - 1) set $2 = _temp2 delete _temp1 delete _temp2 ecleq 4 #ecleq lambda beta ra dec Conversion code from ecliptic to #equatorial coordinates. ra in hours, all others #in degrees. Results are approximate. set _lambda = $1*pi/180 set _beta = $2*pi/180 set _sinb = sin(_beta)*cos(0.3927) set _sinb = _sinb - cos(_beta)*sin(_lambda)*sin(0.3927) set _sinb = ((_sinb > 1) ? 1 : _sinb) set _sinb = ((_sinb < -1) ? -1 : _sinb) set _b = asin(_sinb) set _cosl = cos(_beta)*cos(_lambda)/cos(_b) set _sinl = (cos(_beta)*sin(_lambda)*cos(0.3927)) set _sinl = (_sinl + sin(_beta)*sin(0.3927))/cos(_b) set _cosl = ((_cosl > 1) ? 1 : _cosl) set _cosl = ((_cosl < -1) ? -1 : _cosl) set _l = acos(_cosl) set _l = ((_sinl > 0) ? _l : (2.*pi - _l)) set $4 = _b*180./pi set _l = _l*180./pi set $3 = ((_l > 360) ? (_l - 360) : _l) set $3 = $3/15. delete _sinb delete _sinl delete _cosl delete _l delete _b delete _lambda delete _beta eqgal 4 #eqgal ra dec, l, b. Conversion code. ra in hours, all others #in degrees set _alpha = $1*15*pi/180. set _delta = $2*pi/180. set _sinb = sin(_delta)*cos(1.0926) set _sinb = _sinb - cos(_delta)*sin(_alpha - 4.9262)*sin(1.0926) set _sinb = ((_sinb > 1) ? 1 : _sinb) set _sinb = ((_sinb < -1) ? -1 : _sinb) set _b = asin(_sinb) set _cosl = cos(_delta)*cos(_alpha - 4.9262)/cos(_b) set _sinl = (cos(_delta)*sin(_alpha - 4.9262)*cos(1.0926)) set _sinl = (_sinl + sin(_delta)*sin(1.0926))/cos(_b) set _cosl = ((_cosl > 1) ? 1 : _cosl) set _cosl = ((_cosl < -1) ? -1 : _cosl) set _l = acos(_cosl) set _l = ((_sinl > 0) ? _l : (2.*pi - _l)) set $4 = _b*180./pi set _l = _l*180./pi + 33. set $3 = ((_l > 360) ? (_l - 360) : _l) delete _sinb delete _sinl delete _cosl delete _l delete _b delete _alpha delete _delta galeq 4 #eqgal l b ra dec Conversion code. ra in hours, all others #in degrees set _l = $1*pi/180. set _b = $2*pi/180. set _sind = sin(_b)*cos(1.0926) set _sind = _sind + cos(_b)*sin(_l - 0.5760)*sin(1.0926) set _sind = ((_sind > 1) ? 1 : _sind) set _sind = ((_sind < -1) ? -1 : _sind) set _dec = asin(_sind) set _sinra = (cos(_b)*sin(_l - 0.5760)*cos(1.0926) - \ sin(_b)*sin(1.0926))/cos(_dec) set _cosra = cos(_b)*cos(_l - 0.5760)/cos(_dec) set _cosra = ((_cosra > 1) ? 1 : _cosra) set _cosra = ((_cosra < -1) ? -1 : _cosra) set _ra = acos(_cosra) set _ra = ((_sinra > 0) ? _ra : (2.*pi - _ra)) set $4 = _dec*180./pi set _ra = _ra*180./pi + 282.25 set _ra = ((_ra > 360) ? (_ra - 360) : _ra) set $3 = _ra/15. delete _sind delete _sinra delete _cosra delete _l delete _b delete _ra delete _dec galgrid 0 #Puts up equal area grid #This command only makes sense in landscape mode location 2000 31000 6500 25918 ltype 1 limits -180 180 -90 90 do _i = -75,75,15 { relocate (-180.*cosd($_i)) $_i draw (180.*cosd($_i)) $_i } relocate -182 0 putlabel 4 0 relocate 0 0 putlabel 5 180 relocate 182 0 label 360 do _i = -180, 180, 30 { set _temp = -90, 90 set _line = $_i*cosd(_temp) connect _line _temp } delete _temp delete _line ltype 0 galgridhalf 1 #Puts up equal area grid in upper or lower half #This command only makes sense in portrait mode if ($1 == 1) {location 2000 31000 19161 29988} if ($1 == 2) {location 2000 31000 2778 13605} ltype 1 limits -180 180 -90 90 do _i = -75,75,15 { relocate (-180.*cosd($_i)) $_i draw (180.*cosd($_i)) $_i } relocate -182 0 putlabel 4 0 relocate 0 0 putlabel 5 180 relocate 182 0 label 360 do _i = -180, 180, 30 { set _temp = -90, 90 set _line = $_i*cosd(_temp) connect _line _temp } delete _temp delete _line ltype 0 galpoints 2 #Puts points on equal area grid; galpoints l b set _gl = ($1 - 180.)*cosd($2) points _gl $2 delete _gl galplot 0 #Emulates galplot command readcoor ra dec galgrid eqgal ra dec l b ltype 0 galpoints l b expand 0.8 relocate 80 -70 label $(dimen(l)) points plotted echo $(dimen(l)) points plotted expand 1.2 relocate -144 108 aitoff 4 #converts l and b to x y in aitoff coordinates: #aitoff l b x y set _rl = 90. - 0.5*$1 set _den = sqrt(1 + cosd($2)*cosd(_rl)) set $3 = 2.*cosd($2)*sind(_rl)/_den set $4 = 90.*sind($2)/_den set $3 = 180. - 90.*$3 delete _rl delete _den invaitoff 4 #converts x y in aitoff coordinates to l, b #invaitoff x y l b set _xp = 1. - $1/180. set _yp = $2/90. set _xp = _xp*_xp set _yp = _yp*_yp set _sin2b = _yp*(2 - _xp - _yp) set _b = sqrt(_sin2b) set _b = (_b >= 1)?(0.9999):(_b) set _b = (_b <= -1)?(-0.9999):(_b) set _b = asind(_b) set $4 = ($2 > 0)?(_b):(-_b) set _h = (1 - _xp - _yp)/cosd($4) set _h = (_h >= 1)?(0.9999):(_h) set _h = (_h <= -1)?(-0.9999):(_h) set _l = acosd(_h) set _l = ($1 < 180)?(_l):(-_l) set $3 = 180 - 2*_l foreach i (_xp _yp _sin2b _b _h _l) {delete $i} aitoff1 4 #converts l and b to x y in reverse aitoff coordinates: #aitoff l b x y set _rl = 90. - 0.5*$1 set _den = sqrt(1 + cosd($2)*cosd(_rl)) set $3 = 2.*cosd($2)*sind(_rl)/_den set $4 = 90.*sind($2)/_den set $3 = 180. - 90.*$3 set $3 = 360. - $3 delete _rl delete _den aitoffgrid 0 #Puts up equal area grid in Aitoff #This command only makes sense in landscape mode location 2000 31000 6500 25918 ltype 1 limits 0 360 -90 90 do _i = -75,75,15 { set _l = 0, 360 set _b = _l set _b = $_i aitoff _l _b _x _y connect _x _y } relocate -2 0 putlabel 4 0 relocate 180 0 putlabel 5 180 relocate 362 0 label 360 do _i = 0, 360, 30 { set _b = -90, 90 set _l = _b set _l = $_i aitoff _l _b _x _y connect _x _y } delete _l delete _b delete _x delete _y ltype 0 aitoffgridhalf 1 #Puts up equal area grid in upper or lower half #This command only makes sense in square mode if ($1 == 1) {location 3000 30000 17809 31309} if ($1 == 2) {location 3000 30000 1425 14925} ltype 1 ltype 1 limits 0 360 -90 90 do _i = -75,75,15 { set _l = 0, 360 set _b = _l set _b = $_i aitoff _l _b _x _y connect _x _y } relocate -2 0 putlabel 4 360 relocate 180 0 putlabel 5 180 relocate 362 0 label 0 do _i = 0, 360, 30 { set _b = -90, 90 set _l = _b set _l = $_i aitoff _l _b _x _y connect _x _y } delete _l delete _b delete _x delete _y ltype 0 aitoffgrid3 1 #Puts up equal area grid in three panels. #This command only makes sense in square mode if ($1 == 1) {location 7500 25500 1425 10425} if ($1 == 2) {location 7500 25500 11867 20867} if ($1 == 3) {location 7500 25500 22309 31309} ltype 1 ltype 1 limits 0 360 -90 90 do _i = -75,75,15 { set _l = 0, 360 set _b = _l set _b = $_i aitoff _l _b _x _y connect _x _y } relocate -2 0 putlabel 4 360 relocate 180 0 putlabel 5 180 relocate 362 0 label 0 do _i = 0, 360, 30 { set _b = -90, 90 set _l = _b set _l = $_i aitoff _l _b _x _y connect _x _y } delete _l delete _b delete _x delete _y ltype 0 aitoffgrid1 0 #Puts up equal area grid in Aitoff. Reverse sense #This command only makes sense in landscape mode location 2000 31000 6500 25918 ltype 1 limits 0 360 -90 90 do _i = -75,75,15 { set _l = 0, 360 set _b = _l set _b = $_i aitoff _l _b _x _y connect _x _y } relocate -2 0 putlabel 4 360 relocate 180 0 putlabel 5 180 relocate 362 0 label 0 do _i = 0, 360, 30 { set _b = -90, 90 set _l = _b set _l = $_i aitoff _l _b _x _y connect _x _y } delete _l delete _b delete _x delete _y ltype 0 aitoffplot 0 #Like galplot, but with aitoff projection readcoor ra dec aitoffgrid eqgal ra dec l b ltype 0 aitoff l b x y points x y expand 0.8 relocate 260 -88 label $(dimen(l)) points plotted echo $(dimen(l)) points plotted expand 1.2 relocate 36 108 aitoffplot1 3 #Like galplot, but with aitoff projection #file is in zcore output format. #aitoffplot1 file1 rmin rmax data $1 read {l 1 b 2 r 4} aitoffgrid ltype 0 aitoff l b x y set x = x if (r > $2 && r < $3) set y = y if (r > $2 && r < $3) points x y expand 0.8 relocate 260 -88 label $(dimen(x)) points plotted echo $(dimen(x)) points plotted expand 1.2 relocate 36 108 cum_norm 2 #Obtain normalized cumulant of $1 after sorting; #place in $2 sort < $1 > set _temp1 = $1 cumulate _temp1 _temp2 set $2 = _temp2/sum($1) set _temp2 = (($2 < 1)?($2):(1.)) set $2 = ((_temp2 > 0)?(_temp2):(0.)) delete _temp1 delete _temp2 log_fac 1 # Use Stirling's formula to calculate a log factorial. set $0 = $1 + 1 set $0=$1 == 0 ? 0 : \ 0.5*ln(2*pi) -$0 + ($0-0.5)*ln($0) + ln (1 + 1/(12*$0) + \ 1/(288*$0**2) - 139/(51840*$0**3) - 571/(2488320*$0**4)) chi_dist 3 #Returns, in $3, chi^2 distribution of vector $1 with $2 dof define _x ($2/2. - 1.) define _gamma (log_fac($_x)) set $3 = exp( -0.5*$1 + ($2/2. - 1)*ln(0.5*$1) - $_gamma)/2. define _x delete define _gamma delete emacs_hist ## edit the history list using an external editor. # if you have an environment variable EDITOR it'll be used, # if you don't, we'll use emacs define EDITOR : if(!$?EDITOR) { DEFINE EDITOR emacs } del1 MACRO all 0 100000 # define "all" from buffer !rm -f .all~ MACRO WRITE all .all~ WRITE STANDARD Editing History Buffer\n !sed -e 's/^all[ ]*/ /' -e 's/^..//' .all~ > .all !$EDITOR .all DELETE HISTORY 0 100000 # empty history buffer RESTORE .all !rm -f .all edt 0 #Brings up emacs editor for sm emacs_hist include 1 #like input, but the is put into the history buffer; all #else is cleaned out. del1 DELETE HISTORY 0 10000 restore $1 nint 1 #Compute nearest integer set _temp = abs($1) set _temp = (_temp - int(_temp) <=0.5) ? int(_temp) : int(_temp + 1) set $0 = ($1 > 0 )? _temp : -_temp define _temp delete maxmin 3 #Compute maximum and minimum value of vector $1, put in 2 #and 3 set _temp = $1 sort {_temp} define _size ( dimen(_temp) - 1 ) define $2 ( _temp[$_size] ) define $3 ( _temp[0] ) define _size delete cosd 1 #Returns cosine of the argument, in degrees set $0 = cos(pi*$1/180.) sind 1 #Returns sine of the argument, in degrees set $0 = sin(pi*$1/180.) tand 1 #Returns tangent of the argument, in degrees set $0 = tan(pi*$1/180.) acosd 1 #Returns arc-cosine of the argument, in degrees set $0 = 180.*acos($1)/pi asind 1 #Returns arc-sine of the argument, in degrees set $0 = 180.*asin($1)/pi atand 1 #Returns arc-tangent of the argument, in degrees set $0 = 180.*atan($1)/pi atan2d 2 #Returns arc-tangent of the argument, in degrees set $0 = 180.*atan2($1,$2)/pi sinh 1 #Returns hyperbolic sine set $0 = (e**$1 - e**(-$1))/2. cosh 1 #Returns hyperbolic cosine set $0 = (e**$1 + e**(-$1))/2. tanh 1 #Returns hyperbolic tangent set $0 = (e**$1 - e**(-$1))/(e**$1 + e**(-$1)) asinh 1 #Returns inverse hyperbolic sine. set $0 = ln($1 + sqrt($1**2 + 1)) ave 1 #Return average of argument set $0 = (sum($1)/dimen($1)) sd 1 #Return standard deviation of argument set $0 = (sqrt(sum(($1 - (sum($1)/dimen($1)))**2)/dimen($1))) rms 1 #Return root mean square of argument set $0 = (sqrt(sum($1**2)/dimen($1))) rmsclip 3 #rmsclip x n y Return rms of x in y, clipping at n sigma. set _temp = $1 define _size ($(dimen(_temp))) do i = 1, 10 { define _sigma ($(rms(_temp))) define _mean ($(ave(_temp))) set _temp = _temp if \ (_temp < $_mean + $2*$_sigma && _temp > $_mean - $2*$_sigma) define $3 $_sigma if ($(dimen(_temp)) == $_size) {return} define _size ($(dimen(_temp))) } sgn 1 #Return sign of argument set $0 = ($1 < 0) ? -1 : 1 signx 1 #Return sign of argument, for a scalar if ($1 > 0 ) {set $0 = 1} if ($1 < 0 ) {set $0 = -1} if ($1 == 0 ) {set $0 = 0} integral 3 #Integrate $2 over $1 (not evenly spaced); return $3 #$1 should be sorted define _sum 0 define _temp (dimen($2) - 1) do _i = 1, $_temp { define _ip ($_i - 1) define _sum ($_sum + 0.5*($2[$_ip] + $2[$_i])*($1[$_i] - $1[$_ip])) } define _ip delete define _temp delete define $3 $_sum define _sum delete vfield4 4 #vfield4 x y vx vy Vector field at points x, y, #arrows are defined by vx, vy (in units of x axis) set _amp = sqrt($3*$3 + $4*$4) set _angle = (_amp >0)?($3/_amp):0 set _angle = acosd(_angle) set _angle = (($4 > 0)?(_angle):(360. - _angle)) vfield2 $1 $2 _amp _angle define _amp delete define _angle delete greyscale1 14 # Draw a grey-scale image. Usage: greyscale1 [npx npy maxweight dmargin] # where npx and npy are the number of points in x and y, # and maxweight is the largest lweight to use. dmargin is used # to fiddle how close the points get to the edge of the box # The greylevels to use are set by glevels # This is a variant on the greyscale macro of Robert; points outside # the glevels range are not plotted. if(!is_set(glevel_vec,3)) { echo Please set levels with the `glevels' command return } if(dimen(glevel_vec) < 1) { echo You need at least one level return } define NX image if(!$?NX) { echo Image is not defined return } if(!$?1){ define 1 137 } if(!$?2){ define 2 137 } if(!$?3){ define 3 10 } if(!$?4){ define 4 0.25 } set _npx=$1 set _npy=$2 set _maxweight=$3 define lweight | set _pxy=0,_npx*_npy-1 set _py=int((_pxy+0.000001)/_npx) set _px=_pxy-_npx*_py define X0 image define X1 image define Y0 image define Y1 image if($X0<$fx1){ define X0 $fx1 } if($X1>$fx2){ define X1 $fx2 } if($Y0<$fy1){ define Y0 $fy1 } if($Y1>$fy2){ define Y1 $fy2 } set _px=($X1 - $X0)/(_npx + 2*$4 - 0.5)*\ ((_px+$4) concat (_px+0.5+$4)) + $X0 set _py=($Y1 - $Y0)/(_npy + 2*$4 - 0.5)*\ ((_py+$4) concat (_py+0.5+$4)) + $Y0 set _pz=image(_px,_py) ptype 0 0 minmax _Min _Max set _lev=glevel_vec sort{_lev} #if(_lev[0]<$_Min){ # set _lev=$_Min concat _lev concat $_Max #}else{ # set _lev=_lev concat $_Max #} # set _lev=_lev if(_lev>=$_Min && _lev<=$_Max) if(dimen(_lev)>2){ set _weight=0.51,_maxweight,(_maxweight-0.51)/(dimen(_lev)-1) }else{ lweight $3 points _px _py if($(_lev[1]) > _pz && _pz >= $(_lev[0])) } # lweights < 0.5 are special if (dimen(_lev) > 2) { do i=0,dimen(_lev) - 2 { points _px _py if($(_lev[$i+1]) > _pz && _pz >= $(_lev[$i])) }} lweight $lweight foreach 0 {_px _py _pxy _pz _npx _npy _maxweight _lev _weight}{ delete $0 } foreach 0 {X0 X1 Y0 Y1 _Max _Min lweight}{ define $0 delete } sgtrans 4 #Galactic to supergalactic; sgtrans l b sgl sgb set _sgx = {-0.7357425689697266 0.6772612929344177 0} set _sgy = {-.074553773384391064 -0.080991465669853824 0.9939225912094116} set _sgz = {0.6731452991992128 0.7312711606134599 0.1100812554359436} set _nx = cosd($2)*cosd($1) set _ny = cosd($2)*sind($1) set _nz = sind($2) set _x = _nx*$(_sgx[0]) +_ny*$(_sgx[1]) +_nz*$(_sgx[2]) set _y = _nx*$(_sgy[0]) +_ny*$(_sgy[1]) +_nz*$(_sgy[2]) set _z = _nx*$(_sgz[0]) +_ny*$(_sgz[1]) +_nz*$(_sgz[2]) set _sgb = asind(_z) set _p = _x/cosd(_sgb) set _p = (_p >= 1)?(0.9999):(_p) set _p = (_p <= -1)?(-0.9999):(_p) set _sgl = acosd(_p) set _sgl = (_y >= 0)?(_sgl):(360. - _sgl) set $3 = _sgl set $4 = _sgb foreach 0 {_sgx _sgy _sgz _nx _ny _nz _p _sgb _sgl} { define $0 delete} gstrans 4 #Supergalactic to Galactic; gstrans sgl sgb l b set _sgx = {-0.7357425689697266 0.6772612929344177 0} set _sgy = {-.074553773384391064 -0.080991465669853824 0.9939225912094116} set _sgz = {0.6731452991992128 0.7312711606134599 0.1100812554359436} set _nx = cosd($2)*cosd($1) set _ny = cosd($2)*sind($1) set _nz = sind($2) set _x = _nx*$(_sgx[0]) +_ny*$(_sgy[0]) +_nz*$(_sgz[0]) set _y = _nx*$(_sgx[1]) +_ny*$(_sgy[1]) +_nz*$(_sgz[1]) set _z = _nx*$(_sgx[2]) +_ny*$(_sgy[2]) +_nz*$(_sgz[2]) set _sgb = asind(_z) set _p = _x/cosd(_sgb) set _p = (_p >= 1)?(0.9999):(_p) set _p = (_p <= -1)?(-0.9999):(_p) set _sgl = acosd(_p) set _sgl = (_y >= 0)?(_sgl):(360. - _sgl) set $3 = _sgl set $4 = _sgb foreach 0 {_sgx _sgy _sgz _nx _ny _nz _p _sgb _sgl} { define $0 delete} j1 1 #The first spherical Bessel function set $0 = ($1 > 0.01)?((sin($1) - $1*cos($1))/$1/$1):($1/3.) hprint 2 # Print a help file on topic $1 to file $2 if(!$?help) { define help : set _help = '$help' define help DELETE } else { set _help = '$help' } !/bin/rm -f $2 _hprint $1 $2 DELETE _hprint _hprint 2 # work routine for hprint do 3=1,1000 { set _i = index(_help,' ') if(_i < 0) { set _i = 0 } define _dir (substr(_help,0,_i)) ! test -f $_dir""$1 if($exit_status == 0) { # found it echo Saving history entry for $1 to $2 !cp $_dir""$1 $2 return } set _help = substr(_help,_i + 1,0) if(_i == 0) { # no more directories on path echo Can't find an entry for $1 return } } calc_lb 6 #Calculates r, l and b given x, y, and z #calc_lb x y z r l b set $4 = sqrt($1**2 + $2**2 + $3**2) set _sinb = ($4 > 0)?($3/$4):(0) set _sinb = (_sinb < 1)?(_sinb):(0.999999) set _sinb = (_sinb > -1)?(_sinb):(-0.999999) set $6 = asind(_sinb) set _cosl = ($4 > 0 && ($6 != 90 && $6 != 270))?($1/$4/cosd($6)):(0) set _cosl = (_cosl < 1)?(_cosl):(0.999999) set _cosl = (_cosl > -1)?(_cosl):(-0.999999) set $5 = acosd(_cosl) set $5 = ($2 > 0)?($5):(360. - $5) delete _sinb delete _cosl calc_xyz 6 #Calculates x, y, z given r, l, b #calc_xyz r l b x y z set $4 = $1*cosd($2)*cosd($3) set $5 = $1*sind($2)*cosd($3) set $6 = $1*sind($3) pielabel 7 #Puts on axes for a redshift pieplot #Arguments are: #1 Coordinate system (e/g/s = equatorial, galactic, supergalactic) #2 (1/2) Arc coordinate is first or second #3-6 lower and upper limits in first and second coordinate #7 Maximum redshift (km/s) #Input is assumed to be km/s, and will be labelled as such, #unless variable kms is defined and set to zero. LOCATION 3500 31000 3500 31000 define cut $1 define direction $2 define llim1 $3 define ulim1 $4 define llim2 $5 define ulim2 $6 if ($llim1 > $ulim1) { if ('$cut' == 'e') {define llim1 ($llim1 - 24)} else { define llim1 ($llim1 - 360)}} define vmax $7 expand 1.0001 angle 0 if ('$cut' == 'e') {#A cut in Equatorial coordinates define llim1 ($llim1*15) define ulim1 ($ulim1*15)} if ($direction == 1) { #The arc dimension is the first variable define meanl (0.5*($llim1 + $ulim1)) define meanb (0.5*(cosd($llim2) + cosd($ulim2))) define meanl ($meanl*$meanb) define maxtheta ($ulim1*$meanb - $meanl) } if ($direction == 2) { #The arc dimension is the second variable define meanb (0.5*($llim2 + $ulim2)) define maxtheta ($ulim2 - $meanb) } define xmax ($vmax*sind($maxtheta)) define ymax ($vmax*cosd($maxtheta)) limits $(-$vmax/2.) $($vmax/2.) 0 $vmax location $gx1 $gx2 $gy1 $gy2 if ($xmax > 0.5*$vmax) { #We need to reset limits so #as to make things fit location $gx1 $($gx2-1000) $gy1 $($gy2-1000) limits $(-$xmax) $xmax 0 $(2.*$xmax)} if ($?twopanel) { if ($twopanel) {location $gx1 $gx2 $gy1 $($gy1 + 0.67*($gy2-$gy1))}} relocate 0 0 draw $xmax $ymax define ticksize (lg($vmax/10)) if ($ticksize - int($ticksize) > 0.3) { define ticksize (10.**(int($ticksize) + 1)) } else { define ticksize (2*(10.**(int($ticksize)))) } define int ($vmax/60) do tick = $ticksize, $vmax, $ticksize { define x ($tick*sind($maxtheta)) define y ($tick*cosd($maxtheta)) relocate $x $y define x1 ($x - $int*cosd($maxtheta)) define y1 ($y + $int*sind($maxtheta)) draw $x1 $y1 define x1 ($x + 2.*$int) relocate $x1 $y putlabel 6 $tick } relocate 0 0 define nxmax (-$xmax) draw $nxmax $ymax do tick = $ticksize, $vmax, $ticksize { define x (-$tick*sind($maxtheta)) define y ($tick*cosd($maxtheta)) relocate $x $y define x1 ($x + $int*cosd($maxtheta)) define y1 ($y + $int*sind($maxtheta)) draw $x1 $y1 } set theta = -$maxtheta, $maxtheta, 0.1 set arc_x = $vmax*sind(theta) set arc_y = $vmax*cosd(theta) connect arc_x arc_y define tick_interval 10 if ('$cut' == 'e' && $direction == 1) { define tick_interval 15} if ($direction == 1) { define start_tick ($tick_interval*(int(($llim1 - 0.001)/$tick_interval) + 1)) do tick = $start_tick, $ulim1, $tick_interval { define theta ($tick*$meanb - $meanl) if ('$cut' == 'e') {define theta $(-$theta)} define x1 ($vmax*sind($theta)) define y1 ($vmax*cosd($theta)) relocate $x1 $y1 define x1 (0.98*$vmax*sind($theta)) define y1 (0.98*$vmax*cosd($theta)) draw $x1 $y1 define x1 (1.04*$vmax*sind($theta)) define y1 (1.04*$vmax*cosd($theta)) relocate $x1 $y1 if ('$cut' == 'e') { define ra ($tick/15.) if ($ra < 0) {define ra (24 + $ra)} putlabel 5 $ra^{\it h} } else { if ($tick < 0) {define tick (360 + $tick)} putlabel 5 $tick^o }}} else { define start_tick ($tick_interval*(int(($llim2 - 0.001)/$tick_interval) + 1)) do tick = $start_tick, $ulim2, $tick_interval { define theta ($tick - $meanb) define x1 ($vmax*sind($theta)) define y1 ($vmax*cosd($theta)) relocate $x1 $y1 define x1 (0.98*$vmax*sind($theta)) define y1 (0.98*$vmax*cosd($theta)) draw $x1 $y1 define x1 (1.04*$vmax*sind($theta)) define y1 (1.04*$vmax*cosd($theta)) relocate $x1 $y1 putlabel 5 $tick^o}} define midx (0.5*$vmax*sind($maxtheta) + 0.1*$vmax*cosd($maxtheta)) define midy (0.5*$vmax*cosd($maxtheta) - 0.1*$vmax*sind($maxtheta)) if ($xmax > 0.5*$vmax) { define midx (0.85*$xmax*sind($maxtheta) + 0.25*$xmax*cosd($maxtheta)) define midy (0.75*$xmax*cosd($maxtheta) - 0.15*$xmax*sind($maxtheta))} relocate $midx $midy angle $(90 - $maxtheta) if (!$?kms) { putlabel 5 Redshift (km s^{-1}) } else { if ($kms == 0) { putlabel 5 Redshift } else {putlabel 5 Redshift (km s^{-1})}} define labelpoint ($vmax*1.1) relocate 0 $labelpoint angle 0 if ('$cut' == 'e' && $direction == 1) { putlabel 5 Right Ascension \alpha} if ('$cut' == 'e' && $direction == 2) { putlabel 5 Declination \delta} if ('$cut' == 'g' && $direction == 1) { putlabel 5 Galactic Longitude {\it l}} if ('$cut' == 'g' && $direction == 2) { putlabel 5 Galactic Latitude {\it b}} if ('$cut' == 's' && $direction == 1) { putlabel 5 Supergalactic Longitude {\it L}} if ('$cut' == 's' && $direction == 2) { putlabel 5 Supergalactic Latitude {\it B}} #define midx (-0.5*$vmax*cosd($maxtheta) - 0.03$vmax) #define midy (0.5*$vmax*sind($maxtheta)) #relocate $midx $midy define bottom (-0.1*$vmax) relocate 0 $bottom #Switch back for labelling if ($llim1 < 0) { if ('$cut' == 'e') {define llim1 ($llim1 + 24)} else { define llim1 ($llim1 + 360)}} if ('$cut' == 'e' && $direction == 1) { putlabel 5 $llim2^o < \delta < $ulim2^o} if ('$cut' == 'e' && $direction == 2) { putlabel 5 $($llim1/15)^{\it h} < \alpha < $($$ulim1/15)^{\it h}} if ('$cut' == 'g' && $direction == 1) { putlabel 5 $llim2^o < {\it b} < $ulim2^o} if ('$cut' == 'g' && $direction == 2) { putlabel 5 $llim1^o < {\it l} < $ulim1^o} if ('$cut' == 's' && $direction == 1) { putlabel 5 $llim2^o < {\it B} < $ulim2^o} if ('$cut' == 's' && $direction == 2) { putlabel 5 $llim1^o < {\it L} < $ulim1^o} pieplot 3 #Plots redshift points; requires pielabel to have #been run. pieplot l b v, where l and b are in #Galactic coordinates. foreach vec (ra dec) {set $vec local} set _l = $1 set _b = $2 set _v = $3 if ('$cut' == 's') { #A cut in Supergalactic coordinates sgtrans _l _b sgl sgb set _l = sgl set _b = sgb} if ('$cut' == 'e') {#A cut in Galactic coordinates galeq _l _b ra dec set _l = ra*15 set _b = dec} set ltemp = _l set btemp = _b foreach x {_l _b _v} { if ($llim1 < $ulim1) {set $x = $x if (ltemp > $llim1 && ltemp <= $ulim1 && \ btemp > $llim2 && btemp <=$ulim2 && _v <= $vmax) } else { set $x = $x if ((ltemp > $llim1 || ltemp <= $ulim1) && \ btemp > $llim2 && btemp <=$ulim2 && _v <= $vmax)}} if ($direction == 1) { #The arc dimension is the first variable set theta = (_l*$meanb - $meanl) if ($llim1 > $ulim1) {set theta = (_l > $llim1)?((_l - 360)*$meanb - $meanl):theta} if ('$cut' == 'e') {set theta = -theta} } if ($direction == 2) { #The arc dimension is the second variable set theta = (_b - $meanb) } set x = _v*sind(theta) set y = _v*cosd(theta) points x y pieplot2 4 #Plots redshift points; requires pielabel to have #been run. pieplot l b r v, where l and b are in #Galactic coordinates. Like pieplot, but this #draws lines between distance (r) and redshift (v). set _l = $1 set _b = $2 set _r = $3 set _v = $4 if ('$cut' == 's') { #A cut in Supergalactic coordinates sgtrans _l _b sgl sgb set _l = sgl set _b = sgb} if ('$cut' == 'e') {#A cut in Galactic coordinates galeq _l _b ra dec set _l = ra*15 set _b = dec} set ltemp = _l set btemp = _b set _r = (_r > $vmax && _v < $vmax)?($vmax):_r set _v = (_v > $vmax && _r < $vmax)?($vmax):_v set rtemp = _r set vtemp = _v foreach x {_l _b _v _r} { if ($llim1 < $ulim1) {set $x = $x if (ltemp > $llim1 && ltemp <= $ulim1 && \ btemp > $llim2 && btemp <=$ulim2 && (rtemp <= $vmax || vtemp <= $vmax)) } else { set $x = $x if ((ltemp > $llim1 || ltemp <= $ulim1) && \ btemp > $llim2 && btemp <=$ulim2 && (vtemp <= $vmax || rtemp <= $vmax))}} if ($direction == 1) { #The arc dimension is the first variable set theta = (_l*$meanb - $meanl) if ($llim1 > $ulim1) {set theta = (_l > $llim1)?((_l - 360)*$meanb - $meanl):theta} if ('$cut' == 'e') {set theta = -theta} } if ($direction == 2) { #The arc dimension is the second variable set theta = (_b - $meanb) } set x = _v*sind(theta) set y = _v*cosd(theta) set xr = _r*sind(theta) set yr = _r*cosd(theta) points x y points xr yr do i = 0, (dimen(x) - 1) { relocate (x[$i]) (y[$i]) ltype 0 if (_v[$i] < _r[$i]) {ltype 2} draw (xr[$i]) (yr[$i]) } ltype 0 circlelabel 5 #Puts on axes for a redshift circleplot #Arguments are: #1 Coordinate system (e/g/s = equatorial, galactic, supergalactic) #2 (1/2) angular coordinate is first or second coordinate #3, 4 lower and upper limits in cut #5 Maximum redshift (km/s) define cut $1 define direction $2 define llim $3 define ulim $4 define vmax $5 expand 1.0001 angle 0 if ($direction == 1) { #The arc dimension is the first variable define meanb (0.5*(cosd($llim) + cosd($ulim))) } if ($direction == 2) { #The arc dimension is the second variable if ('$cut' == 'e') { define llim ($llim*15) define ulim ($ulim*15)} define meanb (0.5*($llim + $ulim)) } limits $(-$vmax) $($vmax) $(-$vmax) $vmax location 3500 31000 3500 31000 set theta = 0, 360, 0.1 set x_arc = $vmax*cosd(theta) set y_arc = $vmax*sind(theta) connect x_arc y_arc define tick_interval 30 do tick = 0, $(360 - $tick_interval), $tick_interval { define x1 ($vmax*cosd($tick)) define y1 ($vmax*sind($tick)) relocate $x1 $y1 define x1 (0.98*$vmax*cosd($tick)) define y1 (0.98*$vmax*sind($tick)) draw $x1 $y1 define x1 (1.04*$vmax*cosd($tick)) define y1 (1.04*$vmax*sind($tick)) relocate $x1 $y1 if ('$cut' == 'e' && $direction == 1) { define ra ($tick/15.) putlabel 5 $ra^{\it h} } else { putlabel 5 $tick^o }} relocate 0 $(1.1*$vmax) putlabel 5 Radius of circle $vmax km s^{-1} angle 0 define ticksize (lg($vmax/10)) if ($ticksize - int($ticksize) > 0.3) { define ticksize (10.**(int($ticksize) + 1)*2) } else { define ticksize (10.**(int($ticksize) + 1)) } set theta = 0, 360, 0.1 do tick = $ticksize, $($vmax - $ticksize), $ticksize { set arc_x = $tick*cosd(theta) set arc_y = $tick*sind(theta) ltype 2 lweight 0.5 connect arc_x arc_y relocate 0 $tick expand 0.7 putlabel 5 $tick} expand 1.001 lweight 2 ltype 0 relocate 0 $(-1.1*$vmax) if ('$cut' == 'e' && $direction == 1) { putlabel 5 Right Ascension \alpha, $llim^o < \delta < $ulim^o} if ('$cut' == 'e' && $direction == 2) { putlabel 5 Declination \delta, $llim^o < \alpha < $ulim^o} if ('$cut' == 'g' && $direction == 1) { putlabel 5 Galactic Longitude {\it l}, $llim^o < {\it b} < $ulim^o} if ('$cut' == 'g' && $direction == 2) { putlabel 5 Galactic Latitude {\it b}, $llim^o < {\it l} < $ulim^o} if ('$cut' == 's' && $direction == 1) { putlabel 5 Supergalactic Longitude {\it L}, $llim^o < {\it B} < $ulim^o} if ('$cut' == 's' && $direction == 2) { putlabel 5 Supergalactic Latitude {\it B}, $llim^o < {\it L} < $ulim^o} circleplot 3 #Plots redshift points; requires circlelabel to have #been run. circleplot l b v, where l and b are in #Galactic coordinates. set _l = $1 set _b = $2 set _v = $3 if ('$cut' == 's') { #A cut in Supergalactic coordinates sgtrans _l _b _sgl _sgb set _l = _sgl set _b = _sgb} if ('$cut' == 'e') {#A cut in Galactic coordinates galeq _l _b _t1 _t2 set _l = _t1*15 set _b = _t2} set temp = _l if ($direction == 1) {set temp = _b} foreach x {_l _b _v} { set $x = $x if (temp > $llim && temp <= $ulim \ && _v <= $vmax) } if ($direction == 1) { #The arc dimension is the first variable set theta = (_l*(cosd(_b)/$meanb)) } if ($direction == 2) { #The arc dimension is the second variable set theta = _b } set x = _v*cosd(theta) set y = _v*sind(theta) points x y limits 0 1 0 1 relocate -0.02 1.02 expand 1.001 label $(dimen(x)) Galaxies limits $(-$vmax) $($vmax) $(-$vmax) $vmax foreach i (_l _b _v _sgl _sgb _ra _dec) {delete $i} circle_exclude 1 #After circleplot, draws lines corresponding to #Galactic latitude excluded regions. set _l = 0, 360, 0.1 set _b = _l*0 + $1 galeq _l _b _rap _decp set _ra = _rap if (abs(_decp - $meanb) < 1) #There are two such values: define _ra1 $(_ra[0]*15) relocate 0 0 ltype 0 draw ($vmax*cosd($_ra1)) ($vmax*sind($_ra1)) define _ra2 $(_ra[$(dimen(_ra) - 1)]*15) relocate 0 0 ltype 0 draw ($vmax*cosd($_ra2)) ($vmax*sind($_ra2)) foreach i (_ra _rap _decp) {delete $i} foreach i (_ra1 _ra2) {define $i delete} log 1 #Replaces argument with its log, base 10, and -99 if less than zero. define _verbose $verbose VERBOSE 0 set _temp = $1 if ($1 <= 0) if ($(dimen(_temp)) > 0) { echo Warning: Setting $(dimen(_temp)) values to -99!} set _temp = ($1 > 0) ? (lg($1)) : (-99) set $1 = _temp delete _temp VERBOSE $_verbose delete _verbose antilog 1 #Replaces argument with its antilog. Inverse of log #operation set $1 = 10.**$1 cleanlog 2 #Replaces argument with its log, base 10, and $2 if less than zero. define _verbose $verbose VERBOSE 0 set _temp = $1 if ($1 <= 0) if ($(dimen(_temp)) > 0) { echo Warning: Setting $(dimen(_temp)) values to $2!} set _temp = ($1 > 0) ? (lg($1)) : ($2) set $1 = _temp delete _temp VERBOSE $_verbose delete _verbose dex 1 #Inverse of the log command; replaces argument with it raised to #the power of 10. -99's are replaced with 0. set $1 = ($1 ==-99)?(0):(10**$1) selfunct 16 #Function selfunct r phi alpha beta xstar vlummin #returns the selection function at values of r. if ($?3 == 1) {set _alpha = $3} if ($?4 == 1) {set _beta = $4} if ($?5 == 1) {set _xstar = $5} if ($?6 == 1) {set _vlummin = $6} set _x = ($1 == 0)?(1):($1/3000.) define _xmin (_vlummin/3000.) set _selfunct = (($_xmin/_x)**(2.*_alpha)) * \ ((1 + $_xmin*$_xmin/_xstar/_xstar)**_beta)/\ ((1 + _x*_x/_xstar/_xstar)**_beta) set $2 = ($1 > _vlummin)? (_selfunct) : 1 delete _selfunct delete _x define _xmin delete lumfunct 17 #Function lumfunct r lum alpha beta xstar vlummin n1 #returns the luminosity function at values of #r = sqrt(L/4 pi fmin). if ($?3 == 1) {set _alpha = $3} if ($?4 == 1) {set _beta = $4} if ($?5 == 1) {set _xstar = $5} if ($?6 == 1) {set _vlummin = $6} if ($?7 == 1) {set _n1 = $7} set _x = ($1 == 0)?(1):($1/3000.) define _xmin (_vlummin/3000.) set _selfunct = (($_xmin/_x)**(2.*_alpha)) * \ ((1 + $_xmin*$_xmin/_xstar/_xstar)**_beta)/\ ((1 + _x*_x/_xstar/_xstar)**_beta) set $2 = 2.30*_n1*_selfunct*$1*(_alpha/$1 + _beta*$1/($1*$1 + \ _xstar*_xstar*9e6)) delete _selfunct delete _x define _xmin delete pla 0 ##abbreviation for play del1 play play 0 ##abbreviation for playback del1 playback 0 0 logtick 0 #Most common ticksize for logarithmic ticks ticksize -1 10 -1 10 aitoffdec 1 #Draws a line of constant declination in an aitoff plot. set _ra = 0, 24, 0.01 set _dec = _ra*0 + $1 eqgal _ra _dec l b aitoff l b _x _y findbreak _x _y _xpp _ypp 1 connect _xpp _ypp foreach i (_ra _dec l b _x _y _xpp _ypp) {delete $i} aitoffdec1 1 #Draws a line of constant declination in a reverse aitoff plot. set _ra = 0, 24, 0.01 set _dec = _ra*0 + $1 eqgal _ra _dec l b aitoff1 l b _x _y findbreak _x _y _xpp _ypp 1 connect _xpp _ypp foreach i (_ra _dec l b _x _y _xpp _ypp) {delete $i} findbreak 5 #Sort x y so that xp, yp are continuous define _break 0 set _x = $1 set _y = $2 set _xp = _x set _yp = _y do i = 1, (dimen(_x) - 1) { if (abs(_x[$i] - _x[$i - 1]) > $5) {define _break $i}} if ($_break != 0 ){ do i = 0, (dimen(_x) - $_break - 1) { define j ($i + $_break) set _xp[$i] = _x[$j] set _yp[$i] = _y[$j]} do i = 0, ($_break - 1) { define j (dimen(_x) - $_break + $i) set _xp[$j] = _x[$i] set _yp[$j] = _y[$i]}} set $3 = _xp set $4 = _yp foreach i (_x _y _xp _yp) {delete $i} define _break delete hemidec 1 #Draws a line of constant declination on a pair of #hemispheres (see hemiplot) set ra = 0, 24, 0.01 set dec = ra*0 + $1 eqgal ra dec l b set lp = l if (b > 20) set bp = b if (b > 20) set lp = 360 - lp set lp = lp + $rotate set lp = (lp < 0)? (360 + lp) : lp set lp = (lp > 360)? (-360 + lp) : lp set r = sqrt(2.*(1 - sind(bp))) set x = r*cosd(lp) - sqrt(2) set y = r*sind(lp) findbreak x y xp yp 0.1 connect xp yp set lp = l if (b < -20) set bp = b if (b < -20) set lp = lp - $rotate + 180 set lp = (lp < 0)? (360 + lp) : lp set lp = (lp > 360)? (-360 + lp) : lp set r = sqrt(2.*(1 + sind(bp))) set x = r*cosd(lp) + sqrt(2) set y = r*sind(lp) findbreak x y xp yp 0.1 connect xp yp hemiplot 6 #Sky plot in two hemispheres, equal area #hemiplot l b v rmin rmax rotation_angle set l = $1 set b = $2 set v = $3 define rmin $4 define rmax $5 define rotate $6 set temp = v if (v > $rmin && v <= $rmax) location 2000 31000 6500 25918 define sqrt2 (sqrt(2)) set l = l if (b > 0 && v > $rmin && v < $rmax) set b = b if (b > 0 && v > $rmin && v < $rmax) set l = 360 - l set l = l + $rotate set l = (l < 0)? (360 + l) : l set l = (l > 360)? (-360 + l) : l set r = sqrt(2.*(1 - sind(b))) set x = r*cosd(l) - sqrt(2) set y = r*sind(l) limits -2.8284 2.8284 -1.4142 1.4142 ptype 20 3 expand .3 points x y set theta = 0, 370 set x = sqrt(2)*cosd(theta) - sqrt(2) set y = sqrt(2)*sind(theta) connect x y # set theta = 95,445 set theta = 0, 360 do circle = 20, 80, 20 { set x = sqrt(2*(1 - sind($circle)))*cosd(theta) - sqrt(2) set y = sqrt(2*(1 - sind($circle)))*sind(theta) # define x (-sqrt(2)) # define y (sqrt(2*(1 - sind($circle)))) # relocate $x $y # expand .5 # putlabel 5 $circle # expand 1 ltype 1 connect x y } do theta = 15, 360, 15 { define thetap (360 - $theta) define thetap ($thetap + $rotate) if ($thetap < 0) {define thetap (360 + $thetap)} if ($thetap > 360) {define thetap (-360 + $thetap)} define xp (sqrt(2*1.1)*cosd($thetap) - sqrt(2)) define yp (sqrt(2*1.1)*sind($thetap)) if ($thetap >= 30 && $thetap <=330) { angle (90+$thetap) expand .7 relocate $xp $yp putlabel 5 $theta} if ($theta/30 == int($theta/30) && $theta <= 180) { define xp (sqrt(2*(1 - sind(20)))*cosd($thetap) - sqrt(2)) define yp (sqrt(2*(1 - sind(20)))*sind($thetap)) relocate $xp $yp ltype 1 define thetap ($thetap + 180) define xp (sqrt(2*(1 - sind(20)))*cosd($thetap) - sqrt(2)) define yp (sqrt(2*(1 - sind(20)))*sind($thetap)) draw $xp $yp ltype 0}} expand 1 ltype 0 angle 0 relocate (-sqrt(2)) (sqrt(2)*1.20) putlabel 5 Northern Galactic Hemisphere set l = $1 set b = $2 set v = $3 set l = l if (b < 0 && v > $rmin && v < $rmax) set b = b if (b < 0 && v > $rmin && v < $rmax) set l = l - $rotate + 180 set l = (l < 0)? (360 + l) : l set l = (l > 360)? (-360 + l) : l set r = sqrt(2.*(1 + sind(b))) set x = r*cosd(l) + sqrt(2) set y = r*sind(l) ptype 20 3 expand .3 points x y set theta = 0, 370 set x = sqrt(2)*cosd(theta) + sqrt(2) set y = sqrt(2)*sind(theta) connect x y # set theta = 95,445 set theta = 1, 360 do circle = 20, 80, 20 { set x = sqrt(2*(1 - sind($circle)))*cosd(theta) + sqrt(2) set y = sqrt(2*(1 - sind($circle)))*sind(theta) # define x (sqrt(2)) # define y (sqrt(2*(1 - sind($circle)))) # relocate $x $y # expand .5 # putlabel 5 $circle # expand 1 ltype 1 connect x y } ltype 0 do theta = 15, 360, 15 { define thetap ($theta - $rotate + 180) if ($thetap < 0) {define thetap (360 + $thetap)} if ($thetap > 360) {define thetap (-360 + $thetap)} define xp (sqrt(2*1.1)*cosd($thetap) + sqrt(2)) define yp (sqrt(2*1.1)*sind($thetap)) if ($thetap >= 210 || $thetap <= 150) { relocate $xp $yp angle (90+$thetap) expand .7 putlabel 5 $theta} if ($theta/30 == int($theta/30) && $theta <= 180) { define xp (sqrt(2*(1 - sind(20)))*cosd($thetap) + sqrt(2)) define yp (sqrt(2*(1 - sind(20)))*sind($thetap)) relocate $xp $yp ltype 1 define thetap ($thetap + 180) define xp (sqrt(2*(1 - sind(20)))*cosd($thetap) + sqrt(2)) define yp (sqrt(2*(1 - sind(20)))*sind($thetap)) draw $xp $yp ltype 0}} expand 1 angle 0 relocate (sqrt(2)) (sqrt(2)*1.20) putlabel 5 Southern Galactic Hemisphere relocate 0 -1.4142 putlabel 5 $(dimen(temp)) Galaxies relocate 0 1.5 hemiconvert 5 #Convert from l b to x y for hemisphere plotting. #hemiconvert l b x y rotate set _l = $1 if ($2 > 0) set _b = $2 if ($2 > 0) set _l = 360 - _l set _l = _l + $5 set _l = (_l < 0)? (360 + _l) : _l set _l = (_l > 360)? (-360 + _l) : _l set _r = sqrt(2.*(1 - sind(_b))) set $3 = _r*cosd(_l) - sqrt(2) set $4 = _r*sind(_l) # set _l = $1 if ($2 < 0) set _b = $2 if ($2 < 0) set _l = _l - $5 + 180 set _l = (_l < 0)? (360 + _l) : _l set _l = (_l > 360)? (-360 + _l) : _l set _r = sqrt(2.*(1 + sind(_b))) set _x = _r*cosd(_l) + sqrt(2) set _y = _r*sind(_l) set $3 = $3 concat _x set $4 = $4 concat _y foreach i (_l _b _x _y _r) {delete $i} hard 0 ##abbreviation for hardcopy del1 hardcopy morgage 7 #morgage price down interest term rate rate1 fraction #price is price of house, and down is the down payment #interest is the interest rate in %, and term is the number of years #over which it is paid. rate is morgage payment and #rate1 is total amount paid for house: #payment of morgage + taxes, insurance, maintenance and utilities - #tax deductions #fraction is ratio of total amount paid to morgage company to total morgage set principal = $1 - $2 set price = $1 set interest = $3/100. set term = $4 set $5 = principal*interest/(1 - exp(-interest*term)) set $7 = ($2 + principal*term*interest/(1 - exp(-interest*term)) - $5*term*0.28)/$1 set $6 = $5 + $taxes + price*0.02483 - 0.28*($taxes + $5) set $5=$5/12 set $6=$6/12 lo 0 end pwd 0 del1 !pwd clear 0 del1 !clear default 0 del1 macro read "/u/strauss/sm/default" emacro 0 #edit default macro del1 !emacs -nw /u/strauss/sm/default macro read "/u/strauss/sm/default" sidehist 2 #Draw a histogram on the side of a figure #sidehist x y define _size (dimen($1) - 2) define _startx ($1[0]) define _starty ($2[0]) relocate $_starty $_startx do _i = 1, $_size{ define _nextx ($1[$_i]) define _nexty ($2[$_i]) draw $_nexty $_nextx define _nextx ($1[$($_i + 1)]) draw $_nexty $_nextx } foreach i (_size _startx _starty _i _nextx _nexty) { define $i delete} atof 1 # Convert a string vector to an arithmetic one. Note that # invalid elements (e.g. 'NaN') can be handled by saying # SET NaN = -1e10 SET x = atof(s) # or SET x = (s == 'NaN') ? -1e10 : atof(s) set $0=1,dimen($1) do 2=0,dimen($1)-1 { set $0[$2] = $($1[$2]) } calc_vr 6 #calc_vr l b vx vy vz vr #calculates radial peculiar velocities given #Galactic coordinates and 3-d vpec in SG coordinates. sgtrans $1 $2 _sgl _sgb calc_xyz 1 _sgl _sgb _x _y _z set $6 = $3*_x + $4*_y + $5*_z foreach i (_sgl _sgb _x _y _z) {delete $i} script 1 #Output history file to a script of name 1.sm del1 MACRO all 0 100000 # define "all" from buffer !rm -f .all~ MACRO WRITE all .all~ !sed -e 's/^all[^a-zA-Z]*/ /' -e 's/^..//' .all~ > $1.sm echo Writing out history buffer to $1.sm cdm 6 #Returns CDM power spectrum, with free parameters Gamma, n, and #sigma8, and h #cdm k power gamma n sigma8 h #k is in h Mpc^{-1} #expression from Efsathiou, Frenk, and White define _n $4 set _kp = $1/$3 set $2 = ($1**$_n)/(1 + (6.4*_kp + (3*_kp)**1.5 + \ (1.7*_kp)**2)**1.13)**(2/1.13) #All the work is in setting the normalization. set _logk = -7,7,0.01 set _k = exp(_logk) set _kp = _k/$3 set _power = (_k**$_n)/(1 + (6.4*_kp + (3*_kp)**1.5 + \ (1.7*_kp)**2)**1.13)**(2/1.13) set _y = _k*8/$6 set _integrand = _k*_k*_k*_power*(3*j1(_y)/_y)**2 define _normal ($5*$5*2*pi*pi/sum(_integrand)/0.01) set $2 = $2*$_normal*$6*$6*$6 foreach i (_normal _n) {define $i delete} foreach i (_kp _logk _k _power _y _integrand) {delete $i} eqgal1 6 #Takes ra and dec (hmsdms) and converts to l and b set _ra = $1 + $2/60. + $3/3600. set _dec = abs($4) + abs($5)/60. + abs($6)/3600. set _dec = ($4 < 0 || $5 < 0 || $6 < 0)?(-_dec):(_dec) eqgal _ra _dec _ll _bb print {_ll _bb} foreach i (_ra _dec _ll _bb) {delete $i} wc 1 #Echos to the screen the dimension of the array in question echo The dimension of vector $1 is $(dimen($1)) dectohms 5 #Converts a decimal number to hours, minutes, and seconds. #dechms number sign hours minutes seconds set $2 = sgn($1) set _temp = abs($1) set $3 = int(_temp) set $4 = int(60*(_temp - $3)) set $5 = 3600.*(_temp - $3 - $4/60.) hmstodec 5 #Converts hours, minutes, seconds, to a decimal number. #hmsdec sign hours minutes seconds number set $5 = $1*($2 + $3/60. + $4/3600.) precess 6 # precess ra dec epoch new_epoch ra_new dec_new #This code has not been tested yet. #ra in hours, dec in degrees. define _dt (0.001*($3 - $4)) define _tau (0.001*($3 - 1900)) define _z0 ((0.1117133 + 6.77527E-4*$_tau)*$_dt + \ (1.4656E-4 - 1.31E-6*$_tau)*$_dt*$_dt + 8.726E-5*$_dt**3) define _z ($_z0 + (3.843118E-4 + 3.1998E-6*$_tau)*$_dt*$_dt + \ 1.5514E-6*$_dt**3) define _zj ((0.097189874 - 4.13692E-4*$_tau - \ 1.79381E-6*$_tau*$_tau)*$_dt - \ (2.0687E-4 + 1.79381E-6*$_tau)*$_dt*$_dt - 2.02652E-4*$_dt**3) define _z (90 + 180*$_z/pi) define _z0 (180*$_z0/pi - 90) define _zj (-180*$_zj/pi) set _ra = $1*15 calc_xyz 1. _ra $2 _x _y _z rotate _x _y _z $_z 3 _xp _yp _zp rotate _xp _yp _zp $_zj 1 _x _y _z rotate _x _y _z $_z0 3 _xp _yp _zp calc_lb _xp _yp _zp _r $5 $6 set $5 = $5/15. foreach i (_dt _tau _z0 _z _zj) {define $i delete} foreach i (_ra _x _y _z _xp _yp _zp _r) {delete $i} rotate 8 #rotate x y z theta d xp yp zp #rotate vector x, y, z, by angle theta #around direction d (1, 2, 3), to #find result xp, yp, zp. theta in degrees if ($5 == 1) { set $6 = $1 set $7 = $2*cosd($4) + $3*sind($4) set $8 = -$2*sind($4) + $3*cosd($4)} if ($5 == 2) { set $7 = $2 set $6 = $1*cosd($4) + $3*sind($4) set $8 = -$1*sind($4) + $3*cosd($4)} if ($5 == 3) { set $8= $3 set $6 = $1*cosd($4) + $2*sind($4) set $7 = -$1*sind($4) + $2*cosd($4)} precess1 8 #precess1 hour ramin rasec deg decmin decsec epoch #new_epoch #prints precessed coordinates to the screen #Will have trouble between 0 and -1 degrees #Seems to have a bug. hmstodec 1. $1 $2 $3 _ra set _sign = sgn($4) set _deg = abs($4) hmstodec _sign $_deg $5 $6 _dec precess _ra _dec $7 $8 _newra _newdec dectohms _newra _temp _hour _ramin _rasec dectohms _newdec _sign _deg _decmin _decsec set _deg = _sign*_deg print {_hour _ramin _rasec _deg _decmin _decsec} foreach i (_sign _deg _ra _dec _newra _newdec _hour \ _ramin _rasec _decmin _decsec) {delete $1} make_hist 2 #Plot histogram of $2, binned by $1 histogram $1 (histogram($2:$1)) median_bin1 4 # Create median of $2 in bins containing # $3 elements of $1, put result in $4 set _x = $1 set _y = $2 sort {_x _y} define _nbin (int(dimen(_y)/$3)) set temp = int($3/2.), $(dimen(_y)), $3 set $4 = 1, $_nbin set $4 = 0 do i = 0, $($_nbin - 1) { set $4[$i] = _y[temp[$i]] } foreach i (_x _y _temp) {delete $i} define nbin delete median_bin2 7 # Create median of $2 in bins $3 wide of $1, # put result in $4 (xaxis is $5); interquartile # range in $6, number in bin in $7 set _x = $1 set _y = $2 maxmin _x _max _min define _nbin (int(($_max - $_min)/$3)) set $4 = 1, $_nbin set $4 = 0*$4 set $5 = 0*$4 set $6 = 0*$4 set $7 = 0*$4 do i = 0, $($_nbin - 1) { set _temp = _y if (_x >= $_min + $i*$3 && _x < $_min + ($i + 1)*$3) sort {_temp} set $4[$i] = _temp[$(dimen(_temp)/2.)] set $6[$i] = _temp[0.75*$(dimen(_temp))] - _temp[0.25*$(dimen(_temp))] set $7[$i] = dimen(_temp) set $5[$i] = $_min + ($i + 0.5)*$3 } foreach i (_x _y _temp) {delete $i} foreach i (_nbin _max _min) {define $i delete} sun_to_lg 3 #sun_to_lg l b v #converts heliocentric v to LG-centered, #using Yahil et al 1977. set _temp = $1 - 105. set $3 = $3 + 308.*(cosd(_temp)*cosd($2)*cosd(-7.) + \ sind($2)*sind(-7.)) delete _temp linfit 4 # linear least squares fit # $1 : list of vectors of the linear system # $2 : right side vector # $3 : unknown vector # $4 : variance on $3 # D_$2 : differences between $2 and $2_calculated # the normal equation matrix is in EN (EN_0 EN_1 etc) # # For example, Given Y = a0 + a1 X + a2 U + a3 sin(V) + e, # a is obtained by: # # set sV = sin(v) # set ONE = 0 * X + 1 # set vec = { ONE X U sV } # linfit vec Y a var_a # set e = D_Y # set sig_a = sqrt(var_a) # define n (dimen($1)) define nm (dimen($1)-1) set dimen($3) = $n set dimen($4) = $n do i=0,$nm { set $4[$i] = sum( $($1[$i]) * $2 ) } eqnorm $1 EN qminv EN do i=0,$nm { set $3[$i] = sum( EN_$i * $4 ) } set D_$2 = $2 do i=0,$nm { set D_$2 = D_$2 - $3[$i] * $($1[$i]) } set D2 = sum( D_$2 * D_$2 ) / dimen( D_$2 ) do i=0,$nm { set $4[$i] = D2 * EN_$i[$i] } nwindow 3 #nwindow nx ny n #Puts up a window on an array of nx, ny, in the nth place define _y (abs($2) - int(($3 - 1)/abs($1))) define _x ($3 - abs($1)*int(($3 - 1)/abs($1))) window $1 $2 $_x $_y define _x delete define _y delete nbox 3 #nwindow nx ny n #Puts up a box on an array of nx, ny, in the nth place, with #the labels done right. define _y ($2 - int(($3 - 1)/$1)) define _x ($3 - $1*int(($3 - 1)/$1)) box 0 0 if ($_y == 1) {box 1 0} if ($_x == 1) {box 0 2} define _x delete define _y delete sdsseq 4 #sdsseq lambda eta ra dec #Converts from SDSS survey coordinates to Ra, dec #lambda, eta, dec in degrees, ra in hours #Follows notation of sdss-southern msg 77 define at_surveyCenterRa 185. define at_surveyCenterDec 32.5 define anode ($at_surveyCenterRa - 90.) set _x1 = -sind($1) set _temp = $2 + $at_surveyCenterDec set _y1 = cosd(_temp)*cosd($1) set _z1 = sind(_temp)*cosd($1) calc_lb _x1 _y1 _z1 _r _ra _dec set $3 = (_ra + $anode)/15. set $3 = ($3 > 24)?($3 - 24):($3) set $3 = ($3 < 0)?($3 + 24):($3) set $4 = _dec foreach i (_x1 _y1 _z1 _temp _r _ra _dec) {delete $i} eqsdss 4 #eqsdss ra dec lambda eta #Converts from Ra, dec to SDSS survey coordinates #lambda, eta, dec in degrees, ra in hours #Follows notation of sdss-southern msg 77 define at_surveyCenterRa 185. define at_surveyCenterDec 32.5 define anode ($at_surveyCenterRa - 90.) set _ra = $1*15. set _dec = $2 set _temp = _ra - $anode set _x1 = cosd(_temp)*cosd(_dec) set _y1 = sind(_temp)*cosd(_dec) set _z1 = sind(_dec) calc_lb _y1 _z1 _x1 _r _eta _lambda set $3 = -_lambda set $4 = _eta - $at_surveyCenterDec set $4 = (abs($3) > 90)?($3 + 180):($4) set $3 = ($3 > 90)?(180 - $3):($3) set $3 = ($3 < -90)?(-180 - $3):($3) set $4 = ($4 > 180)?($4 - 360):($4) set $4 = ($4 < -180)?($4 + 360):($4) foreach i (_x1 _y1 _z1 _temp _r _ra _dec _lambda _eta) {delete $i} eqmunu 4 #eqsdss ra dec mu nu #Converts from Ra, dec to SDSS survey coordinates #mu,nu, eta, dec in degrees, ra in hours #Follows notation of sdss-southern msg 77 #see also http://www.astro.princeton.edu/~ivezic/sdssmoc/astrotools.c define at_surveyCenterRa 185. define at_surveyCenterDec 32.5 define anode ($at_surveyCenterRa - 90.) eqsdss $1 $2 lambda eta set eta_stripe = 2.5*(int(eta/2.5 + 0.5)) define ainc (eta_stripe + $at_surveyCenterDec) define ainc (($ainc>180)?($ainc-360):$ainc) echo lambda eta ainc are $(lambda[0]) $(eta_stripe[0]) $ainc set _ra = $1*15. set _dec = $2 set _x1 = cosd(_ra-$anode)*cosd(_dec) set _y1 = sind(_ra-$anode)*cosd(_dec) set _z1 = sind(_dec) set _x2 = _x1 set _y2 = _y1*cosd($ainc) + _z1*sind($ainc) set _z2 = -_y1*sind($ainc) + _z1*cosd($ainc) calc_lb _x2 _y2 _z2 _r _mu _nu set _mu = _mu + $anode set $3 = _mu set $4 = _nu #set $4 = (abs($3) > 90)?($4 + 180):($4) #set $3 = ($3 > 90)?(180 - $3):($3) #set $3 = ($3 < -90)?(-180 - $3):($3) #set $4 = ($4 > 180)?($4 - 360):($4) #set $4 = ($4 < -180)?($4 + 360):($4) foreach i (_x1 _y1 _z1 _x2 _y2 _z2 _r _ra _dec _mu _nu) {delete $i} more 0 MACRO WRITE all junk~ !sed -e 's/^all[^a-zA-Z]*/ /' -e 's/^..//' junk~ > junk !more junk !\rm -f junk* tail 0 MACRO WRITE all junk~ !sed -e 's/^all[^a-zA-Z]*/ /' -e 's/^..//' junk~ > junk !tail junk !\rm -f junk* printif 2 # Print vectors $2... if $1 is true set ll local set ll=($1) define v local foreach v <$2> { set $v local set $v = $v if(ll) } print < $2 > log_poisson 3 #Returns the log_e Poisson distribution of the #arguments. #log_poisson N Nbar answer set $3 = $1*ln($2) - $2 - log_fac($1) aitoffgridlab 2 #Puts labels on aitoffgridlim plot. #Only makes sense if aitoffgridlim has already #been run. #First index indicates whether l should be #labelled as degrees or hours. #Second is whether it runs positive or #negative. #Still with some glitches. do _i = $_bmin, $_bmax, 15 { aitoff $($_lmin - 0.04*$_xrange) $_i _x _y relocate $(_x[0]) $(_y[0]) putlabel 5 $_i } if ($_lmin < $_lmax) { do _i = $_lmin, $_lmax, 30 { define _j $_i if ('$2' == '-') {define _j (360 - $_i)} aitoff $_j $($_bmin - 0.02*$_xrange) _x _y relocate $(_x[0]) $(_y[0]) define _label ($_i - $_offset) if ($_label < 0) {define _label ($_label + 360)} if ('$1' == 'ra' || '$1' == 'RA') { putlabel 5 $($_label/15.) } else { putlabel 5 $_label} }} else { do _i = $_lmin, 330, 30 { define _j $_i if ('$2' == '-') {define _j (360 - $_i)} aitoff $_j $($_bmin - 0.02*$_xrange) _x _y relocate $(_x[0]) $(_y[0]) define _label ($_i - $_offset) if ($_label < 0) {define _label ($_label + 360)} echo $_label if ('$1' == 'ra' || '$1' == 'RA') { putlabel 5 $($_label/15.) } else { putlabel 5 $_label}} do _i = 0, $_lmax, 30 { define _j $_i if ('$2' == '-') {define _j (360 - $_i)} aitoff $_j $($_bmin - 0.02*$_xrange) _x _y relocate $(_x[0]) $(_y[0]) define _label ($_i - $_offset) if ($_label < 0) {define _label ($_label + 360)} echo $_label if ('$1' == 'ra' || '$1' == 'RA') { putlabel 5 $($_label/15.) } else { putlabel 5 $_label}} } delete _x delete _y foreach i (_label) {define $i delete} aitoffgridlim 4 #Puts up equal area grid in Aitoff projection #with limits given by arguments. #square mode is assumed; we assume that this #is appropriate in most cases. location 3500 31000 3500 31000 define _lmin $1 define _lmax $2 define _bmin $3 define _bmax $4 define _wrap 0 if ($_lmin > $_lmax) {define _wrap 1} define _offset (180 - 0.5*($_lmin + $_lmax - 360*$_wrap)) define _lmin ($_lmin + $_offset) define _lmax ($_lmax + $_offset) if ($_lmin > 360) {define _lmin ($_lmin - 360)} if ($_lmax > 360) {define _lmax ($_lmax - 360)} if (abs($_bmin) > abs($_bmax)) { aitoff $_lmin $_bmax _x1 _y1 aitoff $_lmax $_bmax _x2 _y2 } else { aitoff $_lmin $_bmin _x1 _y1 aitoff $_lmax $_bmin _x2 _y2} if ($_bmin < 0 && $_bmax > 0) { aitoff $_lmin 0 _x1 _y1 aitoff $_lmax 0 _x2 _y2} define _xrange (abs(_x1[0] - _x2[0])) define _yrange (abs($_bmin - $_bmax)) if ($_xrange > $_yrange) {define _size $_xrange } else { define _size $_yrange} ltype 2 limits $(180 - $_size/2) $(180 + $_size/2) $_bmin $($_bmin + $_size) do _i = $_bmin, $_bmax, 15 { set _l = $_lmin, $_lmax aitoff _l $_i _x _y connect _x _y } if ($_lmin < $_lmax) { do _i = $_lmin, $_lmax, 30 { set _b = $_bmin, $_bmax aitoff $_i _b _x _y connect _x _y }} else { do _i = $_lmin, 0, 30 { set _b = $_bmin, $_bmax aitoff $_i _b _x _y connect _x _y } do _i = 0, $_lmax, 30 { set _b = $_bmin, $_bmax aitoff $_i _b _x _y connect _x _y }} foreach i (_l _b _x _y _x1 _x2 _y1 _y2) {delete $i} foreach i (_yrange _size) {define $i delete} ltype 0 spectrum_fits 3 # spectrum_fits lambda object # Assumes standard output of spectroscopic pipeline define file_type fits image $1 define NAXIS image define NAXIS1 image define NAXIS2 image if ($NAXIS != 2 || $NAXIS2 != 2) { echo This spectrum is not of the right format; NAXIS = $NAXIS; echo NAXIS2 = $NAXIS2 return} set _i = 0, ($NAXIS1 - 1) set $3 = image[_i,0] #Further checks would be appropriate to see that COEFF's are defined. define COEFF0 image define COEFF1 image define COEFF2 image define COEFF3 image if ($?COEFF0 == 0 || $?COEFF1 == 0) { echo The wavelength solution coefficients are not written to the header. return} if ($?COEFF2 == 0) {define COEFF2 0} if ($?COEFF3 == 0) {define COEFF3 0} define WFITTYPE image if ($?WFITTYPE == 0) { echo I don't know what sort of wavelength coefficients I'm echo looking at; there is no WFITTYPE in the header. return} if ('$WFITTYPE' == 'POLYNOMIAL') { set $2 = $COEFF0 + _i*($COEFF1 + _i*($COEFF2 + _i*$COEFF3))} set newpix = (_i + 1)*2/$NAXIS1 - 1 if ('$WFITTYPE' == 'CHEBYSHEV') { set $2 = $COEFF0 + $COEFF1*cheby1(newpix) + $COEFF2*cheby2(newpix) + \ $COEFF3*cheby3(newpix)} cheby1 1 #Returns Chebyshev Polynomial of order 1 (of type 1) set $0 = $1 cheby2 1 #Returns Chebyshev Polynomial of order 2 (of type 1) set $0 = 2*$1*$1 - 1 cheby3 1 #Returns Chebyshev Polynomial of order 3 (of type 1) set $0 = $1*(4*$1*$1 - 3) specmerge 3 #Read final 1-D spectra output from spectroscopic pipeline #specmerge lambda counts define table_type bintable table "spSpec-51269-0057-1-"$1".fit" read table 'byname' {val[0-1999]} set dimen($3) = 2000 do _i = 0, 1999 { set $3[$_i] = val$_i[0]} set _i = 0, 1999 read table 'byname' {disp.coeff[0-1]} set $2 = 10.**(disp.coeff0[0] + _i*(disp.coeff1[0])) delete _i read_trans 12 # Read a trans file (from a suitable default) for band $1 #This routine is updated for the case of a single asTrans #for the whole set of chips. if(!$?2) { define 2 "$!astrodata/asTrans-$!run.fit" } #All the work is figuring out which chip we're looking at. #This indicates which position ugriz have on the chip. set filter_position = {3 5 1 2 4} table $(5*($col - 1) + filter_position[$1]) "$!2" define INCL image define NODE image local define vars "field a b c d e f dRow0 dRow1 dRow2 \ dRow3 dCol0 dCol1 dCol2 dCol3" define v local foreach v ($vars) { set $v local } read table < $vars > #print { field a b c d e f } #echo table $($1+1) "$!2": $(field) $(a) $(b) $(c) $(d) $(e) $(f) foreach v ($vars) { set trans_$v = $v} # get_mu_nu 4 # Given a field, band, row, and col, create vectors # called $3 $4. E.g. # get_mu_nu 12 2 mu2 nu2 # to set the vectors mu2, nu2 to the r' great-circle positions # of objects in field 12 # Note that rowc, colc have already been defined. DCR # correction is *not* incorporated here. read_trans $2 local set i=do(0,dimen(trans_field)-1) if(trans_field == $1) if(dimen(i) == 0) { user abort "Field $!1 is not present in trans table" } set rowm = rowc$2 + trans_dRow0[i] + trans_dRow1[i]*colc$2 + \ trans_dRow2[i]*(colc$2**2) + trans_dRow3[i]*(colc$2**3) set colm = colc$2 + trans_dCol0[i] + trans_dCol1[i]*colc$2 + \ trans_dCol2[i]*(colc$2**2) + trans_dCol3[i]*(colc$2**3) set $3 = trans_a[i] + trans_b[i]*rowm + trans_c[i]*colm set $4 = trans_d[i] + trans_e[i]*rowm + trans_f[i]*colm # munu_radec 4 #Given mu, nu, calculate ra and dec, both in degrees foreach i (angdiff x2 y2 z2 x1 y1 z1) {set $i local} set angdiff = $1 - $NODE set x2 = cosd(angdiff)*cosd($2) set y2 = sind(angdiff)*cosd($2) set z2 = sind($2) set x1 = x2 set y1 = y2*cosd($INCL) - z2*sind($INCL) set z1 = y2*sind($INCL) + z2*cosd($INCL) set $3 = atan2d(y1, x1) + $NODE set $4 = asind(z1) set $3 = ($3 > 360)?($3 - 360):$3 set $3 = ($3 < 0)?($3 + 360):$3 # radec_munu 4 #Given ra, dec (in degrees), return mu, nu foreach i (angdiff x2 y2 z2 x1 y1 z1) {set $i local} set angdiff = $1 - $NODE set x1 = cosd(angdiff)*cosd($2) set y1 = sind(angdiff)*cosd($2) set z1 = sind($2) set x2 = x1 set y2 = y1*cosd($INCL) + z1*sind($INCL) set z2 =-y1*sind($INCL) + z1*cos($INCL) set $3 = atan2d(y2,x2) + $NODE set $4 = asind(z2) set $3 = ($3 > 360)?($3 - 360):$3 set $3 = ($3 < 0)?($3 + 360):$3 rc_radec 4 #rc_radec frame color ra dec #This routine assumes that rowc, colc exist. #ra, dec in degrees. set mu2 local set nu2 local get_mu_nu $1 $2 mu2 nu2 munu_radec mu2 nu2 $3 $4 # get_rowc_colc 6 # Given a field, band, mu, nu, create vectors # called $5 $6. E.g. # get_rowc_colc 12 2 mu0 nu0 r2 c2 # to set the vectors r2, c2 to the r' row,col positions # of u' objects in field 12 read_trans $2 local set i=do(0,dimen(trans_field)-1) if(trans_field == $1) if(dimen(i) == 0) { user abort "Field $!1 is not present in trans table" } local set det = trans_b[i]*trans_f[i] - trans_c[i]*trans_e[i] local set m = $3 - trans_a[i] local set n = $4 - trans_d[i] set $5 = (m*trans_f[i] - n*trans_c[i])/det set $6 = (-m*trans_e[i] + n*trans_b[i])/det # readline 1 #Reads (the first 40 characters of) an entire line, #and puts the result in string vector $1. read '%[^A]' {$1} return_flag 1 #This routine translates an SDSS flag into English set flag_names = { CANONICAL_CENTER BRIGHT EDGE BLENDED CHILD PEAKCENTER \ NODEBLEND NOPROFILE NOPETRO MANYPETRO NOPETRO_BIG TOO_MANY_PEAKS CR \ MANYR50 MANYR90 BAD_RADIAL INCOMPLETE_PROFILE INTERP SATUR NOTCHECKED \ SUBTRACTED NOSTOKES BADSKY PETROFAINT TOO_LARGE DEBLENDED_AS_PSF \ DEBLEND_PRUNED ELLIPFAINT BINNED1 BINNED2 BINNED4 MOVED DETECTED} do i = 0, (dimen(flag_names)-1) { if (is_set($1,$i)) {print '$(flag_names[$i]) ' {}}} if ($1 != 0) {echo} return_flag2 1 #Translates an SDSS flag2 into English set flag_names = {DEBLENDED_AS_MOVING NODEBLEND_MOVING \ TOO_FEW_DETECTIONS BAD_MOVING_FIT STATIONARY \ PEAKS_TOO_CLOSE MEDIAN_CENTRE LOCAL_EDGE \ BAD_COUNTS_ERROR BAD_MOVING_FIT_CHILD \ DEBLEND_UNASSIGNED_FLUX SATUR_CENTER INTERP_CENTER \ DEBLENDED_AT_EDGE DEBLEND_NOPEAK PSF_FLUX_INTERP \ TOO_FEW_GOOD_DETECTIONS CENTER_OFF_AIMAGE \ DEBLEND_DEGENERATE BRIGHTEST_GALAXY_CHILD \ CANONICAL_BAND AMOMENT_FAINT AMOMENT_SHIFT \ AMOMENT_MAXITER MAYBE_CR MAYBE_EGHOST NOTCHECKED_CENTER} do i = 0, (dimen(flag_names)-1) { if (is_set($1,$i)) {print '$(flag_names[$i]) ' {}}} if ($1 != 0) {echo} return_target 1 #This routine translates an SDSS primTarget flag into #English. set flag_names = {QSO_HIZ QSO_CAP QSO_SKIRT \ QSO_FIRST_CAP QSO_FIRST_SKIRT GALAXY_RED GALAXY \ GALAXY_BIG GALAXY_BRIGHT_CORE ROSAT_A ROSAT_B ROSAT_C \ ROSAT_D STBHB STCARBON STBROWN_DWARF STSUB_DWARF \ STCATY_VAR STRED_DWARF STWHITE_DWARF SERENDIP_BLUE \ SERENDIP_FIRST SERENDIP_RED SERENDIP_DISTANT \ SERENDIP_MANUAL QSO_FAINT GALAXY_REDII ROSAT_E STAR_PN \ QSO_REJECT} do i = 0, (dimen(flag_names) - 1) { if (is_set($1,$i)) {print '$(flag_names[$i]) ' {}}} if ($1 != 0) {echo} read_flag 4 #This routine reads, and returns in English, the flags #associated with an object of a given run number, #column, field, and id. define runn ( sprintf('%06d', $1) ) define datadir "/u/dss/data/$!1/calibChunks/$!2" define rerun 1 define ffield ( sprintf('%04d', $3) ) table 1 "$!datadir/tsObj-$!runn""-$!2""-$!rerun""-$!ffield.fit" read table {id objc_flags objc_flags2 flags[0-4] flags2[0-4]} foreach vec (objc_flags objc_flags2 flags0 flags1 flags2 flags3 \ flags4 flags20 flags21 flags22 flags23 flags24) { set $vec = $vec if (id == $4) } echo objc_flags $(objc_flags[0]) $(objc_flags2[0]) return_flag $(objc_flags[0]) return_flag2 $(objc_flags2[0]) do _i = 0, 4 { echo flags$_i $(flags$_i[0]) $(flags2$_i[0]) return_flag $(flags$_i[0]) return_flag2 $(flags2$_i[0]) } contour_dot 15 #A macro to make contour and dot plots from color-color diagrams. #Arguments are quantities in x and y, and optionally, #binsize in x and y and base contour level. define binsizex 0.02 define binsizey 0.02 if ($?3) {define binsizex $3 define binsizey $4} define mincontour 0.2 if ($?5) {define mincontour $5} define xmin ($fx1) define xmax ($fx2) define ymin ($fy1) define ymax ($fy2) define nbinx (int(abs((($xmax - $xmin)/$binsizex)))) define nbiny (int(abs((($ymax - $ymin)/$binsizey)))) echo $nbinx $nbiny set dimen(temp) = $($nbinx*$nbiny) set temp = temp*0 image ($nbinx,$nbiny) $xmin $xmax $ymin $ymax set image[*,*] = temp set binx = abs(int(($1 - $xmin)/$binsizex)) set biny = abs(int(($2 - $ymin)/$binsizey)) #This should be a one-dimensional quantity unique to #each binx, biny pair. set onedindex = binx*$nbiny + biny set ii = 0, $($nbinx*$nbiny-1) set number_in_bin = histogram(onedindex:ii) #Now we just have to unpack this do i = 0, $($nbinx-1) { #echo $i set j = 0, $($nbiny-1) set image[$i,j] = number_in_bin[$i*$nbiny + j]*$nbinx*$nbiny/dimen(binx) } #The histogram routine puts all objects outside the range in the #very last entry. We can fix this as follows: set image[$nbinx-1,$nbiny-1]=image[$nbinx-2,$nbiny-2] #Can we smooth the contours a bit? set x = 0, $($nbinx -1) set y = 0, $($nbiny -1) do x = 1, $($nbinx-1) { set cut = image[$x,y] smooth cut smoothcut 3 set image[$x,y] = smoothcut } do y = 1, $($nbiny-1) { set cut = image[x,$y] smooth cut smoothcut 3 set image[x,$y] = smoothcut } #Next, we need to assign to each point a contour level. set binx = (binx < 0)?(0):(binx) set binx = (binx > $nbinx -1)?($nbinx-1):(binx) set biny = (biny < 0)?(0):(biny) set biny = (biny > $nbiny -1)?($nbiny-1):(biny) set local_density = image[binx,biny] minmax min max #Draw contours define interval (($max-$min)/9) set levs = ($min+$interval*$mincontour),$max,$interval levels levs contour ptype 0 0 #Dither a bit to get rid of discretization set temp1 = $1 + 0.01*(random(dimen($1)) - 0.5) set temp2 = $2 + 0.01*(random(dimen($2)) - 0.5) set draw_points = (local_density < $min + $interval*$mincontour) points temp1 temp2 if (draw_points) echo The total number of points plotted is $(sum(draw_points)) ctype default foreach i (binx biny onedindex number_in_bin draw_points temp1 temp2) {delete $i} covar 19 #Returns the mean and covariance matrix of up to 9 quantities #Mean is in vector _ave #covariance matrix is in _covar #Its inverse is in _invcov #Square root of the diagonal elements of the covariance matrix is #in _err #The normalized covariance matrix is in _normcov define narg 1 while {$?$narg && $narg < 10} {define narg ($narg + 1)} define narg ($narg - 1) echo We're working with $narg dimensions. #First calculate the mean. set dimen(_ave) = $narg do _i = 0, ($narg-1) { set _ave[$_i] = ave($$($_i+1)) } mdeclare _covar $narg do _i = 0, ($narg-1) { define iarg ($_i+1) do _j = $_i, ($narg-1) { define jarg ($_j+1) set _covar_$_i[$_j] = ave(($$iarg - $(_ave[$_i]))*($$jarg - $(_ave[$_j]))) set _covar_$_j[$_i] = _covar_$_i[$_j] } } set dimen(_err) = $narg do _i = 0, ($narg-1) { set _err[$_i] = sqrt($(_covar_$_i[$_i])) } mdeclare _normcov $narg do _i = 0, ($narg-1) { do _j = $_i, ($narg-1) { set _normcov_$_i[$_j] = _covar_$_i[$_j]/(_err[$_i]*_err[$_j]) set _normcov_$_j[$_i] = _normcov_$_i[$_j] } } mset _invcov _covar minv _invcov trim 1 #Trims a string variable of trailing spaces define end_of_string (INDEX('$1',' ')) if ($end_of_string>0) {define $1 (substr('$1',0,$end_of_string))} return_spflag 1 #Translates a spectral mask bit into English set spflag_names = {NOPLUG BADTRACE BADFLAT BADARC \ MANYBADCOLUMNS MANYREJECTED LARGESHIFT \ ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' \ NEARBADPIXEL LOWFLAT FULLREJECT PARTIALREJECT \ SCATTEREDLIGHT CROSSTALK NOSKY BRIGHTSKY NODATA \ COMBINEREJ} do i = 0, $(dimen(spflag_names) - 1){ if (is_set($1,$i)) {print ' $(spflag_names[$i]) ' {}} } print '\n' {} spectro 0 include "/u/strauss/spectro.sm" play include "/peyton/scr/guarana0/lei/AGNanalysis/spectro1.sm" play xterm_b define idl 1 comoving 1 #The comoving distance in Mpc corresponding to redshift $1. #Multiply by (1+z) for luminosity distance (in a flat universe) #Assumes parameters Omega_m, Omega_l, and h (defaulted to #0.35, 0.65, 0.65 respectively) if (!$?Omega_m) {define Omega_m 0.35} if (!$?Omega_l) {define Omega_l 0.65} if (!$?h) {define h 0.65} define step ($1/5000.) define step (($step>0.01)?(0.01):$step) set _z = 0, $1, $step set integrand = 1/sqrt($Omega_m*(1+_z)**3 + (1 - $Omega_m - $Omega_l)*(1+_z)**2 + $Omega_l) set $0 = (sum(integrand)*$step) set $0 = (3000*$0/$h) lookback 1 #The lookback time (in Megayears) corresponding to redshift $1. #Assumes parameters Omega_m, Omega_l, and h (defaulted to #0.35, 0.65, 0.65 respectively) if (!$?Omega_m) {define Omega_m 0.35} if (!$?Omega_l) {define Omega_l 0.65} if (!$?h) {define h 0.65} define amax (1/(1+$1)) define step ($amax/5000.) define step (($step>0.001)?(0.001):$step) set _a = $step, $amax, $step set _z = 1/_a - 1 set integrand = 1/sqrt($Omega_m*(1+_z)**3 + (1 - $Omega_m - $Omega_l)*(1+_z)**2 + $Omega_l) set integrand = integrand/_a set $0 = (sum(integrand)*$step) set $0 = (1e4*$0/$h) linfit_err 5 # linear least squares fit, with error weighting # $1 : list of vectors of the linear system # $2 : right side vector # $3 : errors on right side vector # $4 : unknown vector (output) # $5 : variance on $3 (output) # D_$2 : differences between $2 and $2_calculated # the normal equation matrix is in EN (EN_0 EN_1 etc) # # For example, Given Y = a0 + a1 X + a2 U + a3 sin(V) + e, # a is obtained by: # # set sV = sin(v) # set ONE = 0 * X + 1 # set vec = { ONE X U sV } # linfit_err vec Y Y_err a var_a # set e = D_Y # set sig_a = sqrt(var_a) # define n (dimen($1)) define nm (dimen($1)-1) set dimen($4) = $n set dimen($5) = $n do i=0,$nm { set $5[$i] = sum( $($1[$i]) * $2/ $3/ $3) } eqnorm_err $1 $3 EN qminv EN do i=0,$nm { set $4[$i] = sum( EN_$i * $5 ) } set D_$2 = $2 do i=0,$nm { set D_$2 = D_$2 - $4[$i] * $($1[$i]) } set D2 = sum( D_$2 * D_$2 ) / dimen( D_$2 ) do i=0,$nm { set $5[$i] = D2 * EN_$i[$i] } eqnorm_err 3 # get the normal equations # for a linear least square problem, with error weighting # $1 : list of vectors # $2 : error array # $3 : matrix (no underscore) of normal equations local define _n (dimen($1)) local define _nm ($_n-1) define i local do i=0,$_nm { set dimen($3_$i) = $_n do j=0,$_nm { set $3_$i[$j] = sum($($1[$i])*$($1[$j])/$2 / $2) } } quit 0 ## quit; don't check DEFINE 2 $verbose VERBOSE 0 del1 DEFINE 1 0 # default value QUIT VERBOSE $2 q 0 ## quit; do check DEFINE 2 $verbose VERBOSE 0 del1 DEFINE 1 0 # default value DEFINE 1 ? {Are you sure? Enter 1 or y to really quit} IF('$1' == '1' || '$1' == 'y') { QUIT } QUIT VERBOSE $2 exit 0 ## quit; don't check DEFINE 2 $verbose VERBOSE 0 del1 DEFINE 1 0 # default value QUIT VERBOSE $2 logout 0 ## quit; don't check DEFINE 2 $verbose VERBOSE 0 del1 DEFINE 1 0 # default value QUIT VERBOSE $2 lo 0 ## quit; don't check DEFINE 2 $verbose VERBOSE 0 del1 DEFINE 1 0 # default value QUIT VERBOSE $2