PRO mstarnight,imgset,thresh

;This program is designed to follow the movement of bright stars over a given 
;time and to observe their relative brightness over this time.
;Modified June 15-17, 1994 by T A Oliynyk.
;	Can now handle more than 50 stars.
;       Function POCALOC added, so that this procedure will work for camera
;       positions at EUREKA, RABBIT LAKE, and CALGARY.
;Modified June 20, 1994 by T A Oliynyk
;	Variable THRESH added. Only stars of magnitude less than THRESH will 
;       be used.   
;Modified June 21, 1994 by T A Oliynyk
; 	Postscript file will be saved in landscape form.
;       Flat field correction applied to the images.

COMMON ffstuff,ffimg

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'
        filename=''
	READ,filename
	psopen,filename,/landscape 
ENDIF
;get times at which images were taken
night = '/jasper/cnsr3_data1/'+imgset
fl = findfile(night+'*',count=nfl)
t=LONARR(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)

	IF i EQ 0 THEN BEGIN
	  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)
  	  exp=(kh.misc.exp_scale*kh.exp.exposure)/1000.	
       	ENDIF	
	t(i)=kh.misc.tm.(0)*3600.+kh.misc.tm.(1)*60.+kh.misc.tm.(2)
ENDFOR

;correct for flat field image
getffimg,yr,mo,da

;get generic dark image to subtract.
rdmeandk,FIX(cam),FIX(exp),dk

index=WHERE((t-MIN(t)) LE 13*3600.,nfl);limit the time span of the images
t=t(index)                             ;to 13 hours.
fl=fl(index)
loc=pocaloc(yr,mo,da)
site=pocasite(yr,mo,da)
t=t(WHERE(t GT 0)) ;shorten time array to correspond to number of files used  

;find positions of stars at the time of each image

hms=STRARR(nfl)

FOR i=0,nfl-1 DO BEGIN
	stephem2,yr,mo,da,t(i),loc(0),loc(1),n,desig,az1,za1,mag,vis,thresh=thresh ;pro to find
				          ;azimuth and zenith angle of stars
	IF i EQ 0 THEN BEGIN      
	  stars=N_ELEMENTS(az1) 
	  az=DBLARR(stars,nfl)
	  za=az
	  rise=REPLICATE(max(t),stars) ;times stars become visible
	  set=REPLICATE(min(t),stars)  ;times stars are not visible
	ENDIF

	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


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
	year=kh.misc.dt.year
	month=FIX(kh.misc.dt.month)
	day=FIX(kh.misc.dt.day)

; procedure to find row and col positions from azimuth and z.ang. 
	position,year,month,day,cam,filt,az(*,i),za(*,i),cl,rw,site=site

	IF i EQ 0 THEN BEGIN 
	   cl1=FLTARR(N_ELEMENTS(cl),nfl)
	   rw1=cl1
	ENDIF
	img=(img-dk)>0
	img=FIX(ROUND(img*(ffimg(filt*256:(filt+1)*256-1,cam*256:(cam+1)*256-1)/100.)))
	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
	print,i+1,nfl,FORMAT='(" Loop #:",I3,"     ","Total # of loops:",I3)'

ENDFOR
;plot only the stars that remain on the CCD for the whole time
ind1=WHERE((cl1 LT 0) OR (cl1 GT 255),n1)
ind2=WHERE((rw1 LT 0) OR (rw1 GT 255),n2)

IF n1 NE 0 THEN BEGIN
	cl1(ind1)=0
ENDIF

IF n2 NE 0 THEN BEGIN
        cl1(ind2)=0
ENDIF

ind=FINDGEN(N_ELEMENTS(cl1(*,0)))

FOR i=0,nfl-1 DO BEGIN
	indd=WHERE(cl1(*,i) EQ 0,count)
	IF count NE 0 THEN ind(indd)=0
ENDFOR

ind=WHERE(ind NE 0)
z=N_ELEMENTS(ind)
difft=FLTARR(z,nfl)
cl1t=difft
rw1t=cl1t

FOR i=0,nfl-1 DO BEGIN
	cl1t(*,i)=cl1(ind,i)
	rw1t(*,i)=rw1(ind,i)
	difft(*,i)=diff(ind,i)
ENDFOR


desig=desig(ind)
cl1=cl1t
rw1=rw1t
diff=difft
num=N_ELEMENTS(diff(*,0))
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,3,2,0,0] ;3 cols, 2 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(30)
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
ENDREP UNTIL(tick(i) EQ FLOAT(roundf))

IF tick(0) EQ 0 THEN BEGIN
	tick=tick(WHERE(tick NE 0))
        tick=tick(1:*)
ENDIF ELSE BEGIN
	tick=tick(WHERE(tick NE 0))
	tick=tick(2:*)
ENDELSE

tin=tick MOD 1

IF tin(0) EQ 0.5 THEN BEGIN
	tickm=tick(WHERE(tin EQ 0.5))
ENDIF ELSE BEGIN
	tickm=tick(WHERE(tin EQ 0.0))
ENDELSE

ticksym=STRTRIM(tickm,2)
ticksym=STRING(ticksym,FORMAT='(f5.1)')

FOR k=0,num-1 DO BEGIN
	;plot of star DN vs. UT (hours)
	PLOT,t,REFORM(dn(k,*))>0,XRAN=[round0*3600.,roundf*3600.],XSTYLE=1$
	,YSTYLE=0,TITLE=fl(0),XTICKV=[tickm*3600],XTICKS=N_ELEMENTS(tickm)-1$
	  , XTICKNAME=[ticksym],XTITLE='UT (hours)'
	;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.9,/DEVICE
	IF !D.NAME EQ 'X' THEN g=GET_KBRD(1)
ENDFOR
IF !D.NAME EQ 'PS' THEN psclose
RETURN
END