deftilde # define tilde TeXstrings macro label \def\tilde\#1{\move0\advance\height{\#1}by50~\#1} 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} 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_lp 1 # postscript file, laser landscape mode # laser_lp filnamee device postlandfile $1 laser_pp 1 # postscript file, laser portrait mode # laser_pp filename device postportfile $1 laser_pp2 1 # postscript file, laser portrait mode - modified size # laser_pp2 filename device postportfile2 $1 laser_sp 1 # postscript file, laser square mode # laser_sp filename device postencap $1 laser_sp2 1 # postscript file, laser square mode - modified bounding box (vertical) # laser_sp2 filename device postencap2 $1 laser_sp3 1 # postscript file, laser square mode - modified bounding box (horizontal) # laser_sp3 filename device postencap3 $1 laser_sp4 1 # postscript file, laser square mode - modified bounding box (horizontal) # laser_sp4 filename device postencap4 $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 11 # put date in lower left hand corner # if argument is given, this will be used instead of sm $date location 0 32768 0 32768 limits 0 32768 0 32768 expand 0.55 relocate 68 68 if ($?1) { putlabel 9 {\ti $1} } else { putlabel 9 {\ti $date} } putfig 1 # put figure label in lower right hand corner # putfig label location 0 32768 0 32768 limits 0 32768 0 32768 expand 0.6 relocate 32700 68 putlabel 7 {\ti $1} 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 random 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} smalldots 1 # smalldots filename # postscript file which makes small dots on SparcPrinter device postscript-smalldots $1 smalldots_l 1 # smalldots_l filename # postland file which makes small dots on SparcPrinter device postland-smalldots $1 smalldots_p 1 # smalldots_p filename # postland file which makes small dots on SparcPrinter device postport-smalldots $1 smalldots_s 1 # smalldots filename # postscript file which makes small dots on SparcPrinter device postscript-smalldots $1 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 } splitwin 11 # splitwin panel # setup for a split window on laser_p if ($?1) { location 8000 25000 25250 30250 } else { location 8000 25000 14500 25000 } 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 stats3 6 # compute mean, variance, skewness, kurtosis given pdf # stats3 x p mean variance skewness kurtosis # p(x)dx=pdf -- integration to one will be enforced define _n (sum($2)) define $3 (sum($2*$1)/$_n) define $4 (sum($2*($1-$$3)**2)/$_n) define $5 (sum($2*($1-$$3)**3)/$_n/$$4**1.5) define $6 (sum($2*($1-$$3)**4)/$_n/$$4**2-3) 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 tail 3 # set one vector equal to the tail of another # tail 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} ullabel 6 # ullabel y pexpand nsides istyle lexpand label # put label and point in upper left corner, at y=y relocate 1.0 $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.0 $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.0 $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 vminmax 13 # find minimum and maximum values of a vector # vminmax vec min [max] set _x=$1 sort {_x} define $2 (_x[0]) if ($?3) {define $3 (_x[dimen(_x)-1])} delete _x xterm ## xterm window, in convenient location device x11 -bg white -geometry 650x525+450+325 xterm_l ## xterm window, same dimensions as landscape postscript device x11 -bg white -geometry 971x725+1200+100 xterm_l2 ## xterm window, same dimensions as landscape postscript device x11 -bg white -geometry 900x450+350+0 xterm_p ## xterm window, same dimensions as portrait postscript device x11 -bg white -geometry 725x971+500+30 xterm_p2 ## xterm window, same dimensions as portrait postscript device x11 -bg white -geometry 550x900+500+30 xterm_s ## square xterm window device x11 -bg white -geometry 800x800+1200+0 xterm_bs ## big square xterm window device x11 -bg white -geometry 950x950+300+0 xterm_ss ## small square xterm window device x11 -bg white -geometry 525x525+525+50 gallat 3 # compute galactic latitudes from ra and dec (P-P region) # gallat ra dec b define decg 27.4 define rag 12.0 + 49.0/60.0 define sing (sin(decg*pi/180.)) define cosg (cos(decg*pi/180.)) set _sinb = $sing*sin($2*pi/180.)+$cosg*cos($2*pi/180.)*cos static int called = 0 ; float decg, rag, sing, cosg, sinb, b ; if (called == 0) { called = 1 ; decg = 27.4 ; rag = 12.0 + 49.0/60.0 ; sing = sin(decg*pi/180.) ; cosg = cos(decg*pi/180.) ; } sinb = sing*sin(dec*pi/180.) + cosg*cos(dec*pi/180.)*cos((ra-rag)*pi/12.0) ; b = 180.*asin(sinb)/pi ; return b ; } 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 lineshade 14 # lineshade x y(x) z(x) [Nfine] # shades the region between curves y and z # Nfine = number of lines drawn between consecutive x-values (default=0) if (!$?4) {define _Nfine 0} else {define _Nfine $4} define _dim dimen($1) set dimen(_xfine)=$(($_dim-1)*$_Nfine+$_dim) define _k 0 do _i = 0,dimen($1)-2,1 { do _j = 0,$_Nfine,1 { set _xfine[$_k] = $1[$_i] + ($_j/($_Nfine+1))*($1[$_i+1]-$1[$_i]) define _k $($_k+1) } } set _xfine[dimen(_xfine)-1] = $1[$_dim-1] set dimen(_yfine)=dimen(_xfine) set dimen(_zfine)=dimen(_xfine) define _k 0 do _i = 0,dimen($1)-2,1 { do _j = 0,$_Nfine,1 { set _yfine[$_k] = $2[$_i] + (_xfine[$_k]-$1[$_i])*($2[$_i+1]-$2[$_i])/($1[$_i+1]-$1[$_i]) set _zfine[$_k] = $3[$_i] + (_xfine[$_k]-$1[$_i])*($3[$_i+1]-$3[$_i])/($1[$_i+1]-$1[$_i]) define _k $($_k+1) } } set _yfine[dimen(_xfine)-1] = $2[$_dim-1] set _zfine[dimen(_xfine)-1] = $3[$_dim-1] do _i=0,dimen(_xfine)-1,1 { relocate $(_xfine[$_i]) $(_yfine[$_i]) draw $(_xfine[$_i]) $(_zfine[$_i]) } delete _xfine delete _yfine define _Nfine delete define _i delete define _j delete define _k delete lineshade2 16 # lineshade2 x y(x) z(x) dx stretch [lineangle] # shades the region between curves y and z. dx = x-spacing between lines. # stretch = (y2-y1)/(x2-x1). lineangle = angle of lines (default=45) if (!$?6) {define _lineangle 45} else {define _lineangle $6} define _dim dimen($1) define _m $(tand($_lineangle)*$5) define _dx $($4/sind($_lineangle)) define _xstart $($1[$_dim-1]) define _xend $($1[0]) do _i = 0,$_dim-1,1 { define _cA $($2[$_i]-$_m*$1[$_i]) define _cB $($3[$_i]-$_m*$1[$_i]) define _xA $(-$_cA/$_m) define _xB $(-$_cB/$_m) if($_xA<$_xstart) {define _xstart $_xA} if($_xB<$_xstart) {define _xstart $_xB} if($_xA>$_xend) {define _xend $_xA} if($_xB>$_xend) {define _xend $_xB} } do _xx = $_xstart,$_xend,$_dx { define _yn0 0 define _ynN 0 define _yn1 0 define _yn2 0 define _c $(-$_m*$_xx) define _ysect $($_m*$1[0]+$_c) if(($_ysect>$2[0] && $_ysect<$3[0]) || ($_ysect>$3[0] && $_ysect<$2[0])) { define _x0 $($1[0]) define _y0 $_ysect define _yn0 1 } define _ysect $($_m*$1[$_dim-1]+$_c) if(($_ysect>$2[$_dim-1] && $_ysect<$3[$_dim-1]) || ($_ysect>$3[$_dim-1] && $_ysect<$2[$_dim-1])) { define _xN $($1[$_dim-1]) define _yN $_ysect define _ynN 1 } do _i = 0,$_dim-2,1 { define _m1 $(($2[$_i+1]-$2[$_i])/($1[$_i+1]-$1[$_i])) define _c1 $($2[$_i] - $_m1*$1[$_i]) if($_m1!=$_m) { define _xsect $(($_c-$_c1)/($_m1-$_m)) } else { define _xsect $($1[$_i]-$_dx) } if($_xsect>=$1[$_i] && $_xsect<$1[$_i+1]) { if($_yn1==0) { define _x1a $_xsect define _y1a $($_m*$_x1a + $_c) define _yn1 1 } else { define _x1b $_xsect define _y1b $($_m*$_x1b + $_c) define _yn1 2 } } } do _i = 0,$_dim-2,1 { define _m2 $(($3[$_i+1]-$3[$_i])/($1[$_i+1]-$1[$_i])) define _c2 $($3[$_i] - $_m2*$1[$_i]) if($_m2!=$_m) { define _xsect $(($_c-$_c2)/($_m2-$_m)) } else { define _xsect $($1[$_i]-$_dx) } if($_xsect>=$1[$_i] && $_xsect<$1[$_i+1]) { if($_yn2==0) { define _x2a $_xsect define _y2a $($_m*$_x2a + $_c) define _yn2 1 } else { define _x2b $_xsect define _y2b $($_m*$_x2b + $_c) define _yn2 2 } } } if($_yn1==0 && $_yn2==0 && $_yn0==1 && $_ynN==1) { define _x1 $_x0 define _y1 $_y0 define _x2 $_xN define _y2 $_yN } if($_yn1==0 && $_yn2==1 && $_yn0==1 && $_ynN==0) { define _x1 $_x0 define _y1 $_y0 define _x2 $_x2a define _y2 $_y2a } if($_yn1==0 && $_yn2==1 && $_yn0==0 && $_ynN==1) { define _x1 $_x2a define _y1 $_y2a define _x2 $_xN define _y2 $_yN } if($_yn1==1 && $_yn2==0 && $_yn0==1 && $_ynN==0) { define _x1 $_x0 define _y1 $_y0 define _x2 $_x1a define _y2 $_y1a } if($_yn1==1 && $_yn2==0 && $_yn0==0 && $_ynN==1) { define _x1 $_x1a define _y1 $_y1a define _x2 $_xN define _y2 $_yN } if($_yn1==1 && $_yn2==1 && $_yn0==0 && $_ynN==0) { define _x1 $_x1a define _y1 $_y1a define _x2 $_x2a define _y2 $_y2a } if($_yn1==2 && $_yn2==0 && $_yn0==0 && $_ynN==0) { define _x1 $_x1a define _y1 $_y1a define _x2 $_x1b define _y2 $_y1b } if($_yn1==0 && $_yn2==2 && $_yn0==0 && $_ynN==0) { define _x1 $_x2a define _y1 $_y2a define _x2 $_x2b define _y2 $_y2b } define _yntot $($_yn0+$_ynN+$_yn1+$_yn2) if($_yntot==2) { relocate $_x1 $_y1 draw $_x2 $_y2 } } define _dx delete define _lineangle delete define _xstart delete define _xend delete define _i delete define _m delete define _c delete define _cA delete define _cB delete define _xA delete define _xB delete define _m1 delete define _c1 delete define _m2 delete define _c2 delete define _xsect delete define _ysect delete define _yn0 delete define _ynN delete define _yn1 delete define _yn2 delete define _yntot delete define _x1a delete define _x1b delete define _y1a delete define _y1b delete define _x2a delete define _x2b delete define _y2a delete define _y2b delete define _x1 delete define _x2 delete define _y1 delete define _y2 delete pointshade 5 # pointshade x y(x) z(x) dx stretch # shades the region between curves y and z. # dx = x-spacing between points. stretch = (y2-y1)/(x2-x1) define _dim dimen($1) define _dx $4 define _dy $($4*$5) do _xx = $($1[0]),$($1[$_dim-1]),$_dx { define _i 0 do _j = 0,$_dim-2,1 { if($_xx>=$1[$_j] && $_xx<$1[$_j+1]) { define _i $_j } } define _yy1 $($2[$_i]+($_xx-$1[$_i])*($2[$_i+1]-$2[$_i])/($1[$_i+1]-$1[$_i])) define _yy2 $($3[$_i]+($_xx-$1[$_i])*($3[$_i+1]-$3[$_i])/($1[$_i+1]-$1[$_i])) if($_yy1>$_yy2) { define _tmp $_yy1 define _yy1 $_yy2 define _yy2 $_tmp } do _yy = $($_yy1+$_dy),$_yy2,$_dy { relocate $_xx $_yy dot } } define _dx delete define _dy delete define _i delete define _j delete define _xx delete define _yy delete define _yy1 delete define _yy2 delete 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} aitoffgrid 11 #Puts up equal area grid in Aitoff #Argument is azimuthal angle at which to break (default 0) #This command only makes sense in landscape mode location 2000 31000 6500 25918 if ($?aitoffsegment) { if ($aitoffsegment == 1) {location 3000 30000 17809 31309} if ($aitoffsegment == 2) {location 3000 30000 1425 14925} } 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 } if (!$?1) {define _break 0} else {define _break $1} #relocate -2 0 #define lower_label (($_break<0)?($_break+360):$_break) #putlabel 4 $lower_label do boo = 0, 360, 60 { relocate $boo -10 if ($?invaitoff) { if ($invaitoff) {relocate $(360-$boo) -10}} define upper_label ($boo+$_break) define upper_label (($upper_label>=360)?($upper_label-360):$upper_label) putlabel 5 $upper_label} #relocate 362 0 #label $lower_label foreach delta (-60 -30 30 60) { aitoff 180 $delta xx yy relocate $(xx[0]+10) $(yy[0]) putlabel 5 $delta} 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 aitoff 15 #converts l and b to x y in aitoff coordinates: #aitoff l b x y if (!$?5) {define _break 0} else {define _break $5} set _rl = ($1-$_break) set _rl = (_rl>360)?(_rl-360):_rl set _rl = (_rl<0)?(_rl+360):_rl set _rl = 90. - 0.5*_rl 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 if ($?invaitoff) { if ($invaitoff) {set $3 = 360. - $3}} delete _rl define _break delete delete _den orthogrid 12 #Puts up orthographic projection grid #orthogrid (default=0,0) #This command only makes sense in square mode if (!$?1) {define _shift_y 0} else {define _shift_y $1} if (!$?2) {define _shift_x 0} else {define _shift_x $2} location 2000 31000 2000 31000 ltype 1 limits -1 1 -1 1 do _i = -75,75,15 { set _l = 0, 360, 1 set _b = _l*0 + $_i orthographic _l _b _xx _yy $_shift_y $_shift_x connect _xx _yy } do _i = 0, 360, 30 { set _b = -90, 90, 1 set _l = _b*0 + $_i orthographic _l _b _xx _yy $_shift_y $_shift_x connect _xx _yy } set _l = 0, 360, 1 set _xx = cosd(_l) set _yy = sind(_l) connect _xx _yy delete _l delete _b delete _xx delete _yy define _shift_x delete define _shift_y delete ltype 0 orthographic 17 #converts l and b to x y in orthographic coordinates: #orthographic l b x y if (!$?5) {define _rotate_y 0} else {define _rotate_y $5} if (!$?6) {define _rotate_x 0} else {define _rotate_x $6} if (!$?7) {define _showback 0} else {define _showback $7} set _rl = ($1-$_rotate_y) set _rl = (_rl>360)?(_rl-360):_rl set _rl = (_rl<0)?(_rl+360):_rl set _x = cosd($2)*cosd(_rl) set _y = sind($2) set _z = cosd($2)*sind(_rl) set _y1 = _y*cosd($_rotate_x) + _z*sind($_rotate_x) set _z1 = _z*cosd($_rotate_x) - _y*sind($_rotate_x) if($_showback==0) { set $3 = _x if(_z1>0) set $4 = _y1 if(_z1>0) } else { set $3 = _x set $4 = _y1 } delete _rl delete _x delete _y delete _z delete _y1 delete _z1 define _rotate_x delete define _rotate_y delete eq2gal 4 #converts equatorial ra and dec to Galactic l and b #eq2gal ra dec l b define eq_incl 62.8717 define eq_meet 282.7667 set _sinb = sind($2)*cosd($eq_incl) - cosd($2)*sind($1-$eq_meet)*sind($eq_incl) set _cosb = sqrt(1-_sinb*_sinb) set _cosl33 = cosd($2)*cosd($1-$eq_meet)/_cosb set _sinl33 = (sind($2)*sind($eq_incl) + cosd($2)*sind($1-$eq_meet)*cosd($eq_incl))/_cosb set $4 = asind(_sinb) set dimen($3)=dimen($1) do i=0,dimen($1)-1,1 { if(_sinl33[$i]>0) { set $3[$i] = acosd(_cosl33[$i]) + 33 } else { set $3[$i] = 360 - acosd(_cosl33[$i]) + 33 } } delete _sinb delete _cosb delete _sinl33 delete _cosl33 define eq_incl delete define eq_meet delete gal2eq 4 #converts Galactic l and b to equatorial ra and dec #gal2eq l b ra dec define eq_incl 62.8717 define eq_meet 282.7667 set _sindec = cosd($2)*sind($1-33)*sind($eq_incl) + sind($2)*cosd($eq_incl) set _cosdec = sqrt(1-_sindec*_sindec) set _cosra282 = cosd($2)*cosd($1-33)/_cosdec set _sinra282 = (cosd($2)*sind($1-33)*cosd($eq_incl) - sind($2)*sind($eq_incl))/_cosdec set $4 = asind(_sindec) set dimen($3)=dimen($1) do i=0,dimen($1)-1,1 { if(_cosra282[$i]<-1) { set _cosra282[$i] = -1 } if(_cosra282[$i]>1) { set _cosra282[$i] = 1 } if(_sinra282[$i]>0) { set $3[$i] = acosd(_cosra282[$i]) + $eq_meet } else { set $3[$i] = 360 - acosd(_cosra282[$i]) + $eq_meet } if($3[$i]>360) {set $3[$i] = $3[$i] - 360} if($3[$i]<0) {set $3[$i] = $3[$i] + 360} } delete _sindec delete _cosdec delete _sinra282 delete _cosra282 define eq_incl delete define eq_meet delete circularaxis 19 #plots a circular axis #circularaxis xscale xcen ycen radius phi1 phi2 stick ltick [phishift] if (!$?9) {define _phishift 0} else {define _phishift $9} q define _lticklength (0.02*$1) define _sticklength (0.01*$1) define _fscale ($4/$1) define _phi1 $5 define _phi2 $6 set _phi = $_phi1+$_phishift,$_phi2+$_phishift,0.1 set _xx = $2 + $4*cosd(_phi) set _yy = $3 + $4*sind(_phi) connect _xx _yy do _theta = -360,720,$7 { if($_theta>=$_phi1 && $_theta<=$_phi2) { define _xx1 ($2 + $4*cosd($_theta+$_phishift)) define _yy1 ($3 + $4*sind($_theta+$_phishift)) define _xx2 ($_xx1-$_sticklength*cosd($_theta+$_phishift)) define _yy2 ($_yy1-$_sticklength*sind($_theta+$_phishift)) relocate $_xx1 $_yy1 draw $_xx2 $_yy2 } } do _theta = -360,720,$8 { if($_theta>=$_phi1 && $_theta<=$_phi2) { define _xx1 ($2 + $4*cosd($_theta+$_phishift)) define _yy1 ($3 + $4*sind($_theta+$_phishift)) define _xx2 ($_xx1-$_lticklength*cosd($_theta+$_phishift)) define _yy2 ($_yy1-$_lticklength*sind($_theta+$_phishift)) relocate $_xx1 $_yy1 draw $_xx2 $_yy2 define _theta1 ($_theta+$_phishift) if($_theta1<0) {define _theta1 ($_theta1+360)} if($_theta1>=360) {define _theta1 ($_theta1-360)} define _theta2 $_theta if($_theta2<0) {define _theta2 ($_theta2+360)} if($_theta2>=360) {define _theta2 ($_theta2-360)} if($_theta1<=180) { if($_theta2<10) { define _xx1 ($2 + $4*cosd($_theta1 + 0.3/$_fscale)) define _yy1 ($3 + $4*sind($_theta1 + 0.3/$_fscale)) } if($_theta2>=10 && $_theta2<100) { define _xx1 ($2 + $4*cosd($_theta1 + 0.6/$_fscale)) define _yy1 ($3 + $4*sind($_theta1 + 0.6/$_fscale)) } if($_theta2>=100) { define _xx1 ($2 + $4*cosd($_theta1 + 1.1/$_fscale)) define _yy1 ($3 + $4*sind($_theta1 + 1.1/$_fscale)) } define _xx2 ($_xx1+$_sticklength*cosd($_theta+$_phishift)) define _yy2 ($_yy1+$_sticklength*sind($_theta+$_phishift)) relocate $_xx2 $_yy2 angle ($_theta1-90) label $_theta2^{\circ} } if($_theta1>180) { if($_theta2<10) { define _xx1 ($2 + $4*cosd($_theta1 - 0.3/$_fscale)) define _yy1 ($3 + $4*sind($_theta1 - 0.3/$_fscale)) } if($_theta2>=10 && $_theta2<100) { define _xx1 ($2 + $4*cosd($_theta1 - 0.6/$_fscale)) define _yy1 ($3 + $4*sind($_theta1 - 0.6/$_fscale)) } if($_theta2>=100) { define _xx1 ($2 + $4*cosd($_theta1 - 1.1/$_fscale)) define _yy1 ($3 + $4*sind($_theta1 - 1.1/$_fscale)) } define _xx2 ($_xx1+1.3*$_lticklength*cosd($_theta+$_phishift)) define _yy2 ($_yy1+1.3*$_lticklength*sind($_theta+$_phishift)) relocate $_xx2 $_yy2 angle ($_theta1+90) label $_theta2^{\circ} } } } delete _phi delete _xx delete _yy define _sticklength define _lticklength define _fscale delete define _phi1 delete define _phi2 delete define _phishift delete define _theta delete define _theta2 delete define _xx1 delete define _yy1 delete define _xx2 delete define _yy2 delete