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} fourwin2 1 # fourwin windownumber # set up for four large square windows on laser_p # like fourwin, but maximizing window size # space at top for labeling if ($1==1) {location 1000 16000 16796 28000} if ($1==2) {location 17000 32000 16796 28000} if ($1==3) {location 1000 16000 4796 16000} if ($1==4) {location 17000 32000 4796 16000} 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 postlandencap $1 laser_pp 1 # postscript file, laser portrait mode # laser_pp filename device postportencap $1 laser_sp 1 # postscript file, laser square mode # laser_sp filename device postencap $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 parabola 4 # find parameters of a parabola fitting 3 points # parabola x y a b -- first 3 elements of (x,y) used # solution is y=a(x-x[0])**2+b*(x-x[0])+y[0] define _x2 ($1[1]-$1[0]) define _x3 ($1[2]-$1[0]) define _y2 ($2[1]-$2[0]) define _y3 ($2[2]-$2[0]) define _denom ($_x3*$_x2**2-$_x3**2*$_x2) define $3 (($_y2*$_x3-$_y3*$_x2)/$_denom) define $4 (($_y3*$_x2**2-$_y2*$_x3**2)/$_denom) foreach var1 {_x2 _x3 _y2 _y3 _denom} {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 threewin 1 # threewin windownumber # set up for three full-page windows on laser_p if ($1==1) {location 9000 25000 22000 30000} if ($1==2) {location 9000 25000 13500 21500} if ($1==3) {location 9000 25000 5000 13000} 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 2.0 $1 ltype 0 relocate 2.3 $1 putlabel 6 $3 ulclinlbl 4 # ullinlbl y ltype ctype label # put label and colored line in upper left corner, at y=y relocate 0.5 $1 ltype $2 ctype $3 draw 2.0 $1 ltype 0 ctype 0 relocate 2.3 $1 putlabel 6 $4 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.0 $1 ltype 0 relocate 7.7 $1 putlabel 4 $3 urclinlbl 4 # urlinlbl y ltype ctype label # put label and colored line in upper right corner, at y=y relocate 9.5 $1 ltype $2 ctype $3 draw 8.0 $1 ltype 0 ctype 0 relocate 7.7 $1 putlabel 4 $4 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 -focus xterm_l ## xterm window, same dimensions as landscape postscript device x11 -bg white -geometry 971x725+300+230 -focus xterm_p ## xterm window, same dimensions as portrait postscript device x11 -bg white -geometry 725x971+500+30 -focus xterm_s ## square xterm window device x11 -bg white -geometry 800x800+325+50 -focus xterm_ss ## small square xterm window device x11 -bg white -geometry 525x525+525+50 -focus 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} 12win2 1 # 12win2 windownumber # set up for 12 square windows on laser_p if ($1==1) {location 2000 11000 24780 31500} if ($1==2) {location 2000 11000 16880 23600} if ($1==3) {location 2000 11000 8980 15700} if ($1==4) {location 2000 11000 1080 7800} if ($1==5) {location 12300 21300 24780 31500} if ($1==6) {location 12300 21300 16880 23600} if ($1==7) {location 12300 21300 8980 15700} if ($1==8) {location 12300 21300 1080 7800} if ($1==9) {location 22600 31600 24780 31500} if ($1==10) {location 22600 31600 16880 23600} if ($1==11) {location 22600 31600 8980 15700} if ($1==12) {location 22600 31600 1080 7800} 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 sun_to_lg 3 #sun_to_lg l b v (from M. Strauss, Oct. 2004) #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 eqgal 4 #eqgal ra dec, l, b. Conversion code. ra in hours, all others #in degrees (from M. Strauss, Oct. 2004, note 1950 RA/DEC) 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 (from M. Strauss, Oct. 2004, note 1950 RA/DEC) 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