Viewing contents of file '../idllib/contrib/buie/fieldobs.pro'
;+
; NAME:
; fieldobs
; PURPOSE:
; Real-time display of standard star locations in the sky
; DESCRIPTION:
;
; CATEGORY:
; Astronomy
; CALLING SEQUENCE:
; fieldobs,file,AMCRIT=amcrit,OBSFILE=obsfile,OBSCODE=obs
; INPUTS:
; file - List of objects (for file format see LOADSTARS) to display
;
; OPTIONAL INPUT PARAMETERS:
;
; KEYWORD INPUT PARAMETERS:
; AMCRIT - Maximum airmass for plotted display (default = 3.0)
; OBSFILE - Location of observatory information file, default is
; /pub/sac0/elgb/data/obscode.dat
; OBSCODE - Observatory code, default = 688, Lowell Observatory
;
; OUTPUTS:
;
; KEYWORD OUTPUT PARAMETERS:
;
; COMMON BLOCKS:
;
; SIDE EFFECTS:
;
; RESTRICTIONS:
;
; PROCEDURE:
;
; MODIFICATION HISTORY:
; 97/05/28, Written by Marc W. Buie, Lowell Observatory
; 98/04/30, MWB, added some items to the plot.
;-
pro fieldobs,file,AMCRIT=amcrit,OBSFILE=obsfile,OBSCODE=obs
if badpar(amcrit,[0,2,3,4,5],0,caller='FIELDOBS: (amcrit) ',default=3.0) then return
if badpar(obsfile,[0,7],0,CALLER='FIELDOBS: (OBSFILE) ', $
DEFAULT='/pub/sac0/elgb/data/obscode.dat') then return
if badpar(obs,[0,2,3],0,caller='FIELDOBS: (obscode) ',default=688) then return
rdobscod,code,alllon,rhosinp,rhocosp,obsname,FILE=obsfile
idx=where(obs eq code,count)
idx=idx[0]
if (count eq 1) then begin
lon = (360.0-alllon[idx])/180.0*!pi
lat = atan(rhocosp[idx],rhosinp[idx])
name=strtrim(obsname[idx],2)
endif else begin
obs = -1
endelse
; Read the special positions file, these won't be in the star catalog
openr,luns,file,/get_lun
sname=''
sinfo=''
sign=''
num=0
while not eof(luns) do begin
readf,luns,sname,rh,rm,rs,sign,dd,dm,ds,ep,sinfo, $
format='(a16,2(1x,i2),1x,f4.1,1x,a1,3(i2,1x),f6.1,1x,a)'
newra = rh + rm/60.0 + rs/3600.0
newdec = dd + dm/60.0 + ds/3600.0
if sign eq '-' then newdec = -1.0*newdec
if num eq 0 then begin
id = sname
ra = newra
dec = newdec
equi = ep
rapm = 0.0
decpm = 0.0
epoch = ep
comment = sinfo
endif else begin
id = [id, sname]
ra = [ra, newra]
dec = [dec, newdec]
equi = [equi, ep]
rapm = [rapm, 0.0]
decpm = [decpm, 0.0]
epoch = [epoch, ep]
comment = [comment,sinfo]
endelse
num=num+1
endwhile
free_lun,luns
ra = ra*15.0/!radeg
dec = dec/!radeg
crital = 0.5*!pi - acos(1.0/amcrit)
while 1 do begin
curtime = systime(1) / 86400.0d0
curjd = curtime + 2440587.5d0
dt = 15.0/60.0/24.0
jd=[curjd,curjd+dt,curjd+dt*2]
moonpos,curjd,moonra,moondec,/radian
sunpos,curjd,sunra,sundec,/radian
mphase,curjd,moonphase
mam = airmass(curjd,moonra,moondec,lat,lon,lst=lst,lha=mlha,alt=malt,/hardie)
sam = airmass(curjd,sunra,sundec,lat,lon,lha=slha,alt=salt,/hardie)
salt=salt*!radeg
malt=malt*!radeg
; Compute time to sunrise
hatojd,!dpi,sunra,lst[0],curjd,jdlclmid
altoha,-12.0/!radeg,sundec,lat,sunhorzha,suntype
jdsset = jdlclmid - (!dpi-sunhorzha)/2.0d0/!dpi
jdsrise = jdlclmid + (!dpi-sunhorzha)/2.0d0/!dpi
jdstr,curjd,0,curtimstr
am=fltarr(3,num)
lha=fltarr(3,num)
for i=0,num-1 do begin
am[*,i] = airmass(jd,ra[i],dec[i],lat,lon,lha=lha0,/hardie)
lha[*,i] = lha0
endfor
z=where(am lt amcrit)
lha = lha*!radeg/15.0
mlha = mlha*!radeg/15.0
slha = slha*!radeg/15.0
; compute solar and lunar elongation angles.
selong = fix(angsep(sunra,sundec,ra,dec) * !radeg + 0.5)
melong = fix(angsep(moonra,moondec,ra,dec) * !radeg + 0.5)
print,' '
print,curtimstr
moondesc=string('moon @ ',malt,moonphase,format='(a,f5.1," deg (",f4.2,")")')
sundesc=string('sun @ ',salt,format='(a,f5.1," deg")')
IF mam[0] ge amcrit THEN BEGIN
xr=max(abs(lha[z]))*[-1,1]
ENDIF ELSE BEGIN
xr=max(abs([lha[z],mlha]))*[-1,1]
ENDELSE
yr=[amcrit,1.0]
plot,[0],[1],/nodata,xr=xr,yr=yr, $
xtitle='W Hour Angle E',ytitle='Airmass',title=curtimstr, $
color=7
setusym,-1
oplot,[mlha],[mam],psym=8,symsize=2.0,color=4
oplot,[slha],[sam],psym=8,symsize=2.0,color=4
setusym,1
oplot,[slha],[sam],psym=8,symsize=0.5,color=4
for i=1.0,amcrit,0.5 do $
oplot,xr,[0,0]+i,color=3,linestyle=1
for i=-3,3,3 do $
oplot,[0,0]+i,yr,color=3,linestyle=1
ss=1.5
if malt gt 0. then mc=9 else mc=2
xyouts,0.02,.02,moondesc,align=0.0,/normal,color=mc,charsize=ss
if salt lt -18.0 then $
xyouts,0.85,.02,sundesc,align=0.5,/normal,color=2,charsize=ss $ ; 2=green
else if salt lt -12.0 then $
xyouts,0.85,.02,'AT '+sundesc,align=0.5,/normal,color=4,charsize=ss $ ; 2=green
else if salt lt -6.0 then $
xyouts,0.85,.02,'NT '+sundesc,align=0.5,/normal,color=4,charsize=ss $ ; 2=green
else if salt lt -0.8 then $
xyouts,0.85,.02,'CT '+sundesc,align=0.5,/normal,color=7,charsize=ss $ ; 2=green
else $
xyouts,0.85,.02,sundesc,align=0.5,/normal,color=9,charsize=ss ; 2=green
if salt lt -12.0 then begin
nightleft = (jdsrise - curjd)*24.0
if nightleft gt 1.0 then color=2 else color=6
nightleft = nightleft*15.0/!radeg
rastr,nightleft,-2,timeleftstr
xyouts,0.85,0.97,timeleftstr+' to N.T.',charsize=ss,color=color,align=0.5,/normal
endif else if curjd lt jdsset then begin
dayleft = (jdsset - curjd)*24.0
if dayleft gt 1.0 then color=9 else color=6
dayleft = dayleft*15.0/!radeg
rastr,dayleft,-2,timeleftstr
xyouts,0.85,0.97,timeleftstr+' to N.T.',charsize=ss,color=color,align=0.5,/normal
endif
print,' # Name HA X Mel Sel'
for i=0,num-1 do begin
if min(am[*,i]) lt amcrit then begin
if lha[0,i] lt 0.0 then side='E' else side='W'
rastr,abs(lha[0,i])*15.0/!radeg,-2,hastr
print,i,id[i],hastr,side,am[0,i],melong[i],selong[i], $
format='(3x,i3,2x,a,3x,a,a,3x,f4.1,3x,i3,2x,i3)'
IF melong[i] lt 30 and malt gt 0.0 THEN pc=1 ELSE pc=9
oplot,lha[*,i],am[*,i],psym=-4,color=pc
if lha[0,i] lt 0.0 then begin
align=0.0
xpos=max(lha[*,i])
ypos=min(am[*,i])
endif else begin
align=1.0
xpos=min(lha[*,i])
ypos=min(am[*,i])
endelse
xyouts,xpos,ypos,' '+strtrim(string(i),2)+' ',align=align,color=pc
endif
endfor
empty ; flush graphics buffer
spawn,'sleep 60'
endwhile
end