PRO starnight,imgset

;This program is designed to follow the movement of bright stars over a given 
;time and to observe their relative brightness over this time.

PRINT,'Do you wish to open a Post Script file? (y/[n])'
post=''
READ,'post'

IF STRLOWCASE(post) EQ 'y' THEN BEGIN
	PRINT,'Name the file, using the complete path'
	READ,filename
	psopen,filename
ENDIF
;get times at which images were taken
RLlat=58.233 & RLlon=256.330
t=LONARR(400)
night = '/jasper/cnsr3_data1/'+imgset
fl = findfile(night+'*',count=nfl)
FOR i=0,nfl-1 DO BEGIN
	;read in headers of images and extract time image was taken 
	rdkihd,fl(i),hb & kh=gethd(hb)
	getdatetime,kh,date,time
	yr=date(0)+1900 &  mo=date(1) & da=date(2) ;same for all images in set
	cam=kh.misc.(0) & filt=kh.misc.(1)
	t(i)=kh.misc.tm.(0)*3600.+kh.misc.tm.(1)*60.+kh.misc.tm.(2)
ENDFOR
t=t(WHERE(t GT 0)) ;shorten time array to correspond to number of files used  
;read in # of stars which will be observed
OPENR,1,'/wesson/user/steele/schee/idl/ybsc4s.dat' ;contains a list of stars and their positions
READF,1,nst ;the first line shows the number of stars in the file
rise=REPLICATE(max(t),nst) ;times stars become visible
set=REPLICATE(min(t),nst) ;times stars become not visible
;rise is initially set to the latest time and set to the earliest time, they
;are changed laterwhen the positions of the stars are calculated
CLOSE,1

;find positions of stars at the time of each image
az=DBLARR(nst,nfl) & za=az
hms=STRARR(400)
FOR i=0,nfl-1 DO BEGIN
	stephem2,yr,mo,da,t(i),RLlat,RLlon,n,desig,az1,za1,mag,vis ;pro to find
				          ;azimuth and zenith angle of stars
	az(*,i)=az1 & za(*,i)=za1
	rise(vis)=(rise(vis) < (t(i))) ;times at which stars are first and last
	set(vis)=(set(vis) > (t(i)))   ;seen 
	print,t(i)
	hms(i)=strsec(t(i))
ENDFOR
;follow only those stars which are visible for at least 1 hour
follow=WHERE((set-rise) GT .75*3600.)
desig=desig(follow) & mag=mag(follow)
siz=N_ELEMENTS(follow)
az=az(follow,*) & za=za(follow,*)
!ORDER=1
diff=FLTARR(siz,nfl) ;array of differences calculated using 'starmag'
tt=SYSTIME(1)
;find the DN of the followed stars with respect to the median BG
img1=INTARR(256,256,nfl)
cl1=FLTARR(50,nfl) & rw1=cl1
FOR i=0,(nfl-1) DO BEGIN ;loop executed once per image; all stars considered
				;in each image
	rdkimg,fl(i),hb,img
	kh=gethd(hb)
	IF kh.exp.(6) EQ 512 THEN img = REBIN(img,256,256,/sample); make 512
							;images into 256
	img1(*,*,i)=img
	PRINT,'Were the images taken at Rabbit lake with the 14 degree offset?'
	offset=''
	READ,offset
	IF STRLOWCASE(offset) EQ 'y' THEN $
	     position,cam,filt,az(*,i),za(*,i),cl,rw,/RL $
	     ELSE position,cam,filt,az(*,i),za(*,i),cl,rw ;procedure to find row and col positions from azimuth and z.ang. 
	cl1(0:N_ELEMENTS(cl)-1,i)=cl/2. & rw1(0:N_ELEMENTS(cl)-1,i)=rw/2.
	cl=cl/2. & rw=rw/2. ;row and col found in 'position' are for 512*512
	starmag,img,cl,rw,d ;procedure to find star DN over median BG
	diff(0:N_ELEMENTS(d)-1,i)=d
ENDFOR
num=N_ELEMENTS(d)
st=N_ELEMENTS(t)
dn=FLTARR(num,nfl)
dn=diff(0:num-1,*)

;plotting results
IF !D.NAME EQ 'X' THEN $ ;otherwise postscript file must be running
	WINDOW,0,XSIZE=1200,YSIZE=900 $
	ELSE DEVICE,XSIZE=24,YSIZE=18
!P.MULTI=[0,4,3,0,0] ;4 cols, 3 rows of graphs per page 
hr0=t(0)/3600 & min0=(t(0)-3600*hr0)/60	; note integer arithmetic!
hrmin0=hr0+min0/60. ; round off first time to nearest 1/2 hour
remain0=FIX(hrmin0+1)-hrmin0 
IF (remain0 GE 0) AND (remain0 LT .25) THEN round0=hr0+1
IF (remain0 GE .25) AND (remain0 LT .75) THEN round0=hr0+.5
IF (remain0 GE .75) THEN round0=hr0
hrf=t(st-1)/3600 & minf=(t(st-1)-3600*hrf)/60; round off last time to 
						;nearest 1/2 hour
hrminf=hrf+minf/60. 
remainf=FIX(hrminf+1)-hrminf 
IF (remainf GE 0) AND (remainf LT .25) THEN roundf=hrf+1
IF (remainf GE .25) AND (remainf LT .75) THEN roundf=hrf+.5
IF (remainf GE .75) THEN roundf=hrf

tick=FLTARR(20)
ticksym=STRARR(20)
i=-1
REPEAT BEGIN ;Define ticks on x-axis - start at first rounded off time (round0)
		;and incriment by 1/2 hr until last time (roundf). 
	i=i+1
	tick(i)=round0+i*.5
	ticksym(i)=STRMID(strsec(3600*tick(i)),0,5) ;convert seconds to HH:MM
ENDREP UNTIL(tick(i) EQ FLOAT(roundf))
gz=WHERE(tick NE 0)
tick=tick(gz) & ticksym=ticksym(gz)
FOR k=0,num-1 DO BEGIN
	;plot of star DN vs. UT (HH:MM)
	PLOT,t,REFORM(dn(k,*))>0,XRAN=[round0*3600.,roundf*3600.],XSTYLE=1$
	,YSTYLE=0,TITLE=fl(0),XTICKV=[tick(1:*)*3600],$
	   XTICKNAME=[ticksym(1:*)],XTITLE='UT (hh:mm)'
	;plot movement of star in image over time in bottom-left of plot
	wlx=!X.CRANGE & wly=!Y.CRANGE ;axis limits - DATA coords.
	; construct lower LH corner coords in DATA coords
	wl0=[wlx(0),wly(0),0.0]
	dvllc=CONVERT_COORD([wlx(0),wly(0),0.0],/DATA,/TO_DEVICE) ; -> DEVICE coords.
	dvulc=CONVERT_COORD([wlx(0),0.7*wly(0)+0.3*wly(1),0.0],/DATA,/TO_DEVICE)
	dviwedge=dvulc(1)-dvllc(1)	; extract subplot edge length - DEVICE coords.
	dvurc=[dvllc(0)+dviwedge,dvllc(1)+dviwedge]	;upper RH corner of sub-plot
	dvspx=REFORM((cl1(k,*)/256.)*(dvurc(0)-dvllc(0))+dvllc(0));coords on plot
	dvspy=REFORM((rw1(k,*)/256.)*(dvurc(1)-dvllc(1))+dvllc(1));of cl & rw
	PLOTS,[dvllc(0),dvurc(0),dvurc(0),dvllc(0)] $
	     ,[dvllc(1),dvllc(1),dvurc(1),dvurc(1)],/DEVICE
	PLOTS,dvspx,dvspy,/DEVICE,PSYM=3
	XYOUTS,dvllc(0)+0.05*dviwedge,dvurc(1)+0.05*dviwedge,STRTRIM(desig(k),2) $
	      ,CHARSIZE=0.4,/DEVICE
	IF !D.NAME EQ 'X' THEN g=GET_KBRD(1)
ENDFOR
IF !D.NAME EQ 'PS' THEN psclose
RETURN
END