PRO dropcals,n,t,files,filtlist,typelist,texplist ; This little routine removes calibration images from the image ; list sequence to simplify cycle identification. imgcateg,files,filtlist,typelist,texplist calfiles=WHERE(typelist EQ 2) one2two,calfiles,files,numrow,camcol keep=setminus(INDGEN(n),numrow) n=N_ELEMENTS(keep) t=t(keep,*) files=files(keep,*) PRINT,n,FORMAT='(I0," image pairs excluding calibration images")' imgcateg,files,filtlist,typelist,texplist RETURN END ;--------------------------------------------------------------------- PRO getcycles,n,t,filtlist,typelist,texplist,cycle,cyclen $ ,debug=debug,cont=cont ; Identify the first images or image pairs of repeated sequences. ; With DOS-controlled image acquisition, typical exposure times ; are 1, 10, and 60 s, with correponding inter-image elapsed times ; of 5, 14, and 64 s; any consistently-occurring elapsed times ; much longer than these probably indicate the end of a data-taking ; sequence and the start of another one. With operator-controlled ; image acquisition, exposure times are similar, but images are ; acquired continuously without dropping back to DOS (so no longer ; intervals between images to pinpoint cycle starts). IF KEYWORD_SET(cont) THEN BEGIN ; Continuous data acquisition requires a different definition of a ; cycle and different methods of identifying the start of a cycle. ; Start with discontinuities in filter sequences. fseq0=filtlist(*,0) ; sequence of Cam 0 filters fseq1=filtlist(*,1) ; sequence of Cam 1 filters ; Remove artificial discontinuities caused by missed images; assume ; only one camera misses an image at a time wneg=WHERE(fseq0 EQ -1,nwneg) IF nwneg GT 0 THEN fseq0(wneg)=fseq1(wneg) wneg=WHERE(fseq1 EQ -1,nwneg) IF nwneg GT 0 THEN fseq1(wneg)=fseq0(wneg) fchng0=fseq0-SHIFT(fseq0,1) ; Cam 0 filter change fchng1=fseq1-SHIFT(fseq1,1) ; Cam 1 filter change ; Now find negative values of fchng# fcycle=WHERE((fchng0 LT 0) AND (fchng1 LT 0),nf) ; ; Now look for changes in image type ; tseq0=typelist(*,0) ; sequence of Cam 0 image types tseq1=typelist(*,1) ; sequence of Cam 1 image types ; As above, remove artificial discontinuities wneg=WHERE(tseq0 EQ -1,nwneg) IF nwneg GT 0 THEN tseq0(wneg)=tseq1(wneg) wneg=WHERE(tseq1 EQ -1,nwneg) IF nwneg GT 0 THEN tseq1(wneg)=tseq0(wneg) tchng0=tseq0-SHIFT(tseq0,1) ; Cam 0 image type change tchng1=tseq1-SHIFT(tseq1,1) ; Cam 1 image type change tcycle=WHERE((tchng0 NE 0) AND (tchng1 NE 0),nt) ; ; Now merge the two definitions to try to catch all possible ; types of cycle. ; CASE 1 OF (nf GT 0) AND (nt GT 0): cycle=[fcycle,tcycle] (nf GT 0) AND (nt EQ 0): cycle=fcycle (nf EQ 0) AND (nt GT 0): cycle=tcycle (nf EQ 0) AND (nt EQ 0): BEGIN MESSAGE,'Can''t identify cycle limits; giving up' $ ,/INFORMATIONAL RETURN END ENDCASE ; STOP cycle=cycle(uniq(cycle,SORT(cycle))) HELP,fcycle,tcycle,cycle ; ENDIF ELSE BEGIN ; dt=(t(*,0)-SHIFT(t(*,0),1))>0 ; ; Try decreasing threshold time below 75 s if no cycles identified. ; thresh=75 REPEAT BEGIN cycle=[WHERE(dt GT thresh),n] ; extra element simplifies finding ; last image in last cycle tcyclen=(cycle-SHIFT(cycle,1))>0 success=(N_ELEMENTS(cycle) GT 2) AND $ (N_ELEMENTS(cycle) GT n/7) AND $ (MAX(tcyclen) LT 7) thresh=thresh-1 ENDREP UNTIL success OR (thresh LT 67) IF success $ THEN MESSAGE,STRING(thresh+1 $ ,FORMAT='("Cycles detected with ",I2,' $ +'"-s threshold")') $ ,/INFORMATIONAL ; If the above technique doesn't work, we probably have continuous ; image acquisition without dropping back to DOS. Try another tack. ; IF (((N_ELEMENTS(cycle) EQ 2) AND (MAX(cycle-[-1,n]) EQ 0)) OR $ ; (N_ELEMENTS(cycle) LT n/7)) THEN BEGIN ; flag=0 ; cycle=tryshift(filtlist,typelist,texplist,flag,debug=debug) ; ENDIF ; IF N_ELEMENTS(flag) EQ 0 THEN flag=0 ; IF flag AND NOT KEYWORD_SET(debug) THEN $ ; PRINT,'Had a spot of trouble sorting out the cycles, you know' ENDELSE IF cycle(0) NE 0 $ ; first cycle SHOULD start at image 0 THEN cycle=[0,cycle] cyclen=(cycle-SHIFT(cycle,1))>0 cycle=cycle(0:N_ELEMENTS(cycle)-2) cyclen=cyclen(1:*) RETURN END ;--------------------------------------------------------------------- PRO opnsaveu,yr,month,day,nwsave,wsave,cwl,rpb,cfsave,rayleighs $ ,bgsubtract,savunit,save00,save01,save02,save03,save04 $ ,save10,save11,save12,save13,save14 @isitdos FOR jsave=0,nwsave-1 DO BEGIN tmp=WHERE(ABS(wsave(jsave)-cwl) LT rpb/2.) tmp=tmp(0) cfsave(tmp)=1 write=Where(StrPos(StrUpCase(saveroot),'D') EQ 0,nwrite) IF nwrite EQ 0 THEN Message,'Writeable save directory not found' fname=saveroot(write)+dd+STRING(yr,month,day,tmp/5,tmp MOD 5 $ ,FORMAT='(3I2.2,2I1,".")') fname=fname(0) ; convert from 1-element array to scalar CASE 1 OF (NOT rayleighs): fname=fname+'cDN' (rayleighs AND (NOT bgsubtract)): fname=fname+'ray' bgsubtract: fname=fname+'bgs' ENDCASE ; (1) IF dos THEN OPENW,tt,fname,/GET_LUN,/BINARY ELSE OPENW,tt,fname,/GET_LUN ; IF dos THEN OPENW,tt,fname,/GET_LUN ELSE OPENW,tt,fname,/GET_LUN savunit(tmp)=tt CASE tmp OF 0: save00=ASSOC(savunit(0),BYTARR(66048)) 1: save01=ASSOC(savunit(1),BYTARR(66048)) 2: save02=ASSOC(savunit(2),BYTARR(66048)) 3: save03=ASSOC(savunit(3),BYTARR(66048)) 4: save04=ASSOC(savunit(4),BYTARR(66048)) 5: save10=ASSOC(savunit(5),BYTARR(66048)) 6: save11=ASSOC(savunit(6),BYTARR(66048)) 7: save12=ASSOC(savunit(7),BYTARR(66048)) 8: save13=ASSOC(savunit(8),BYTARR(66048)) 9: save14=ASSOC(savunit(9),BYTARR(66048)) ENDCASE ; (tmp) ENDFOR ; (jsave) RETURN END ;--------------------------------------------------------------------- ; NEWQL.PRO ; ; Originally written to show images from Polar Camera initial field test ; at Ryning Farm Observatory ; Modified 1993/1/29 by DPS to handle images from image tree structure ; on /jasper/cnsr3_data1 ; Modified 1993/5/11 from SHOW.PRO to display images from both cameras ; in the order they were acquired. ; Modified 1993/6/11 from SEQ.PRO to display large numbers of images ; in compressed form. ; Modified 1993/11/... to show image sequences from Eureka, and to apply ; flat-field corrections and conversion to intensities. ; Modified 1994/03/29 to be compatible with filename formats under DOS. ; PRO newql,year,month,day,hr,min,sec $ ,batch=batch,intens=intens,diffbg=diffbg,debug=debug $ ,save=save,wlsave=wlsave,qlimageplot=qlimageplot COMMON KI_image,h,im COMMON ffstuff,ffimg IF N_PARAMS() LT 6 THEN BEGIN MESSAGE,'Must enter year,month,day,hr,min,sec at least!' RETURN ENDIF setstart=SYSTIME(1) ; Set up OS-specific items @isitdos ; create byte-valued input parameters cam=0B & filt=cam s='' !ORDER=1 ; uncover screen where text is to be printed IF !D.NAME EQ 'X' THEN WDELETE,0 ; If interactive operation desired, prompt user and get input IF NOT KEYWORD_SET(batch) THEN batch=0B fac=1. resp='' IF NOT batch THEN BEGIN PRINT,'Do you want to display images as well as save them? ([Y]/N)' READ,resp dispim=(STRUPCASE(resp) NE 'N') ENDIF ELSE dispim=0B ; (NOT batch) IF dispim THEN BEGIN ; set colour table loadct,13,/silent PRINT,'Enter the scale compression factor to be applied (over normal scale; default 1)' READ,resp IF resp NE '' THEN fac=FIX(resp) ENDIF ; (dispim) ; 1-cycle-per-hour option only works for Eureka data! PCsite=PoCasite(year,month,day) IF dispim AND (PCsite EQ 3) THEN BEGIN PRINT,'View all images (0; default) or only 1 cycle per hour (1)?' READ,resp short=(STRUPCASE(resp) EQ '1') ENDIF ELSE short=0B ; (dispim) IF NOT batch THEN BEGIN PRINT,'Convert DN to intensity? ([Y]/N)' READ,resp rayleighs=(STRUPCASE(resp) NE 'N') ENDIF ELSE BEGIN IF KEYWORD_SET(intens) THEN rayleighs=intens ELSE rayleighs=1B ENDELSE ; (NOT batch) IF rayleighs THEN BEGIN IF NOT batch THEN BEGIN PRINT,'Present backgrounds as differential brightnesses (R/nm)? ([Y]/N)' READ,resp RperA=(STRUPCASE(resp) NE 'N') ENDIF ELSE BEGIN IF KEYWORD_SET(diffbg) THEN RperA=diffbg ELSE RperA=1B ENDELSE ; (NOT batch) ; For now, disable background subtraction (DPS, 940105) ; PRINT,'Do background subtractions? ([Y]/N)' ; READ,resp resp='N' bgsubtract=(STRUPCASE(resp) NE 'N') IF bgsubtract THEN BEGIN PRINT,'Settle for approximate background subtraction? ([Y]/N)' READ,resp approx=(STRUPCASE(resp) NE 'N') ENDIF ; (bgsubtract) ENDIF ELSE BEGIN RperA=0 bgsubtract=0 ENDELSE ; (rayleighs) IF dispim THEN BEGIN PRINT,'Fix the maximum value to be displayed for each camera? (Y/[N])' READ,resp resp=STRUPCASE(resp) ceiling_set=(resp EQ 'Y') IF ceiling_set THEN BEGIN topnum=FLTARR(2) FOR j=0,1 DO BEGIN PRINT,'Enter the ceiling value for camera '+STRTRIM(j,2) READ,top topnum(j)=top ENDFOR ; (j) ENDIF ; (ceiling_set) ENDIF ELSE ceiling_set=0 ; (dispim) IF NOT batch THEN BEGIN PRINT,'Display debugging information? (Y/[N])' READ,resp debug=(STRUPCASE(resp) EQ 'Y') ENDIF ; (NOT batch) IF NOT batch THEN BEGIN IF NOT short THEN BEGIN PRINT,'Save images? ([Y]/N) (This option also creates quick-look data files)' READ,resp imsave=(STRUPCASE(resp) NE 'N') ENDIF ELSE imsave=0 ; (NOT short) ENDIF ELSE IF KEYWORD_SET(save) THEN imsave=save ELSE imsave=0B cfsave=BYTARR(10) savunit=BYTARR(10) IF NOT batch THEN BEGIN IF imsave THEN BEGIN PRINT,'Enter the wavelengths to be saved (nm)' savwav='' READ,savwav nwsave=nwrds(savwav) IF nwsave GT 0 THEN BEGIN wsave=FLTARR(nwsave) FOR jsave=0,nwsave-1 DO wsave(jsave)=FLOAT(getwrd(savwav,jsave)) ENDIF ELSE BEGIN MESSAGE,'No wavelengths entered',/INFORMATIONAL ENDELSE ; (nwsave GT 0) ENDIF ELSE nwsave=0 ; (imsave) ENDIF ELSE BEGIN nwsave=N_ELEMENTS(wlsave) IF nwsave GT 0 THEN BEGIN wsave=wlsave nwsave=N_ELEMENTS(wsave) ENDIF ; (nwsave GT 0) ENDELSE ; (NOT batch) ; make lists of available files of desired type rdilist,2,year,month,day,hr,min,sec,n,t,files ; For production purposes, only run NEWQL if all data from the ; given day are available. IF (hr*60L+min)*60+sec EQ 0 THEN BEGIN mint=MIN(t/3600.,MAX=maxt) IF (mint GT 0.5) OR (maxt LT 23.5) THEN BEGIN MESSAGE,"Full day's data not available; not worth running" $ ,/INFORMATIONAL SETENV,'ACTION=skip' RETURN ENDIF ENDIF IF KEYWORD_SET(qlimageplot) THEN initqlif,year,month,day ; ensure a 4-digit year for display; make reasonable assumptions yr=year IF year LT 1900 THEN BEGIN IF year LT 50 THEN year=year+2000 IF year GT 50 THEN year=year+1900 ENDIF ELSE yr=year MOD 100 ; (year LT 1900) ; create title for image plot IF dispim $ THEN wtitle='Polar Camera '+ymd2date(year,month,day,format='d$ n$ y$') ; define masks for images, specific to dates before/after moonshades getsig,year,month,day,mask0,mask1 sig0=WHERE(mask0) sig1=WHERE(mask1) ; define image portion definitely OUTSIDE imaged area on CCD base=WHERE(SHIFT(dist(256),128,128) GT 150) IF n GT 0 THEN BEGIN ; Reverse the indices of the time and file name arrays as returned ; by rdilist t=TRANSPOSE(t) files=TRANSPOSE(files) ext_list=extension(files) ; Identify and delete calibration files from time and file name arrays, ; so no impossibly long cycles appear dropcals,n,t,files,filtlist,typelist,texplist ENDIF ; (n GT 0) IF n GT 0 THEN BEGIN CASE pocasite(year,month,day) OF 1: cont=1B 2: cont=1B 3: cont=0B ELSE: cont=0B ENDCASE getcycles,n,t,filtlist,typelist,texplist,cycle,cyclen $ ,debug=debug,cont=cont ; If a cycle is too long, investigate. IF NOT batch THEN BEGIN IF MAX(cyclen)/fac GT 5 THEN BEGIN PRINT,'Maximum cycle length is '+STRTRIM(MAX(cyclen),2)+'; check it out!' STOP ENDIF ; (MAX(cyclen) GT 5) ENDIF ; (NOT batch) ; If only hourly cycles are to be shown, reduce cycle and cyclen ; appropriately, to display only dark cycle and the following ; open-shutter cycle each hour. ; ; NOTA BENE: This won't work for continuous data acquisition, such ; as that done at Rabbit lake and Spy Hill. Only good for Eureka! ; IF dispim THEN BEGIN ; IF short THEN BEGIN ; wdkcyc=WHERE(cyclen EQ 2) ; wnextcyc=wdkcyc+1 ; dkcyc=cycle(wdkcyc) ; nextcyc=setminus(cycle(wnextcyc),dkcyc) ; scycle=[dkcyc,nextcyc] ; scyclen=[cyclen(wdkcyc),cyclen(wnextcyc)] ; scyclen=scyclen(0:N_ELEMENTS(scycle)-1) ; ssc=SORT(scycle) ; cycle=scycle(ssc) ; cyclen=scyclen(ssc) ; ENDIF ; (short) ; ENDIF ; (dispim) maxcyclen=MAX(cyclen) cycleparms=STRARR(5,2*maxcyclen); to store image annotation strings fac=1 IF maxcyclen GT 5 THEN fac=2 HELP,cyclen,maxcyclen ; IF dos THEN scrwid=4. ELSE scrwid=5. ; IF dispim AND (maxcyclen GT scrwid*fac) THEN BEGIN ; avoid display problems ; longcyc=STRTRIM(maxcyclen,2) ; with long cycles ; maxninrow=STRTRIM(FIX(scrwid*fac),2) ; minnewfac=STRTRIM(CEIL(maxcyclen/scrwid),2) ; PRINT,'The longest cycle is '+longcyc+' images long; no more than' ; PRINT,maxninrow+' images will fit side by side on the screen.' ; PRINT,'You must change the image scale compression factor.' ; PRINT,'Enter a new factor by which to shrink the images ' $ ; +'(min. '+minnewfac+')' ; READ,fac ; ENDIF ; (dispim AND (maxcyclen GT scrwid*fac)) ; create and label image display window imgsize=256/fac IF dispim THEN BEGIN X_charheight=10 ; try to keep scalings that apply in X windows X_charwidth=8 barwidth=16 ninrow=maxcyclen wwidth=ninrow*imgsize dheight=2*imgsize+barwidth ; run colour bar along bottom display=BYTARR(wwidth,dheight) ; array to hold one cycle of images wheight=dheight+4*X_charheight ENDIF ; (dispim) ; Load the pseudo-flat-field images generated at Eureka, for later ; application to the images before display. getffimg,year,month,day ; Load the Rayleigh-second/DN values for converting to intensities IF rayleighs THEN BEGIN getsens,year,month,day,sens0,sens1,nointeractive=batch sens=[sens0,sens1] ; Load the center wavelengths and rectangular passbands for background ; subtraction rdcwlrpb,cwl,rpb,normunits ENDIF ELSE normunits='DN' ; (rayleighs) ; If images are to be saved, open the associated files to hold them IF nwsave GT 0 THEN BEGIN opnsaveu,yr,month,day,nwsave,wsave,cwl,rpb,cfsave,rayleighs $ ,bgsubtract,savunit,save00,save01,save02,save03,save04 $ ,save10,save11,save12,save13,save14 seqnum=INTARR(10) ENDIF ; (nwsave GT 0) ; The following command is commented out for now. IDL for Windows seems to ; have a limit of 16 files open simultaneously. Rathen than open all the ; tmpNM files at once, they will be opened as needed in 'savqldat.pro'. ; IF imsave THEN iniqlfil,qlunits qlunits=INDGEN(10) ; define dummy values for these variables ; Create the format string for describing background subtractions IF NOT KEYWORD_SET(debug) THEN debug=0B IF debug THEN $ bgform='("Image ",I1,", camera ",I1,";"/" subtract image ",I1' $ +',", camera ",I1)' IF dispim THEN BEGIN IF (!D.NAME EQ 'X') OR (!D.NAME EQ 'WIN') THEN $ WINDOW,0,XSIZE=wwidth,YSIZE=wheight $ ELSE DEVICE,/COLOR,BITS_PER_PIXEL=8 ; PostScript ; keep 'X' char size csz=((1.*X_charheight)/wheight)/((1.*!D.Y_CH_SIZE)/!D.Y_VSIZE) XYOUTS,scl(res=120.,wwidth/2) $ ,scl(res=120.,wheight-3*X_charheight),wtitle $ ,/DEVICE,ALIGNMENT=0.5,CHARSIZE=2.0*csz,WIDTH=twidthn twidthd=CONVERT_COORD([twidthn,0.],/NORMAL,/TO_DEVICE) twidthd=twidthd(0) ; we only need the x-component ; get user's name and display it IF dos THEN username='CNSR3' ELSE SPAWN,'whoami',username XYOUTS,scl(res=120.,wwidth-X_charwidth) $ ,scl(res=120.,wheight-3.9*X_charheight) $ ,'User '+STRUPCASE(username),ALIGNMENT=1.0,/DEVICE,CHARSIZE=csz ; get date/time and display them IF dos THEN BEGIN datetime=SYSTIME() date=STRMID(datetime,4,7)+STRMID(datetime,20,4) ENDIF ELSE BEGIN SPAWN,'date -u',datetime date=STRMID(datetime,4,7)+STRMID(datetime,24,4) ENDELSE time=STRMID(datetime,11,9)+' UT' XYOUTS,scl(res=120.,wwidth-X_charwidth) $ ,scl(res=120.,wheight-2.5*X_charheight) $ ,time,ALIGNMENT=1.0,/DEVICE,CHARSIZE=csz XYOUTS,scl(res=120.,wwidth-X_charwidth) $ ,scl(res=120.,wheight-1.1*X_charheight) $ ,date,ALIGNMENT=1.0,/DEVICE,CHARSIZE=csz ; create and display color bar bar=congrid(BINDGEN(!D.N_COLORS-1,1),ninrow*imgsize,/minus_one) bar=REBIN(bar,ninrow*imgsize,barwidth) help,bar window,2 plot,bar(*,0) k=get_kbrd(1) wdelete,2 display(0,dheight-barwidth)=bar ; mtv,bar,0,0,res=120. ; show image orientation, one time only XYOUTS,scl(3*X_charwidth,res=120.) $ ,scl(barwidth+2*imgsize+3*X_charheight,res=120.) $ ,'N',/DEVICE,ALIGNMENT=0.5,CHARSIZE=csz XYOUTS,scl(res=120.,X_charwidth) $ ,scl(res=120.,barwidth+2*imgsize+1.5*X_charheight),'W' $ ,/DEVICE,ALIGNMENT=0.5,CHARSIZE=csz XYOUTS,scl(res=120.,5*X_charwidth) $ ,scl(res=120.,barwidth+2*imgsize+1.5*X_charheight),'E' $ ,/DEVICE,ALIGNMENT=0.5,CHARSIZE=csz XYOUTS,scl(res=120.,3*X_charwidth) $ ,scl(res=120.,barwidth+2*imgsize),'S',/DEVICE $ ,ALIGNMENT=0.5,CHARSIZE=csz ENDIF ; (dispim) ; Decide where Polar Camera was located PCsite=PoCasite(year,month,day) ; Calculate the azimuth and zenith angle images for each camera col=BINDGEN(256)#REPLICATE(1B,1,256) row=TRANSPOSE(col) cr2az,year,month,day,0,0,col,row,az0,za0,/deg,site=PCsite cr2az,year,month,day,1,0,col,row,az1,za1,/deg,site=PCsite ; Identify image pixels within 2 degrees of the zenith zen0=WHERE(za0 LT 2.) zen1=WHERE(za1 LT 2.) IF bgsubtract THEN BEGIN IF dispim THEN BEGIN WINDOW,3,XSIZE=256,YSIZE=512,XPOS=0,YPOS=0 WSET,0 ENDIF ; (dispim) ; Restrict future attention to image points where data are ; significant - use flat field correction images as guides ; sig0.EU is LONARR(57649) ; sig1.EU is LONARR(58383) IF debug THEN PRINT,'Identifying image interior...' OPENR,s0unit,droot+dd+'sig0.EU',/GET_LUN s0=0L READU,s0unit,s0 sig0=LONARR(s0) READU,s0unit,sig0 CLOSE,s0unit FREE_LUN,s0unit OPENR,s1unit,droot+dd+'sig1.EU',/GET_LUN s1=0L READU,s1unit,s1 sig1=LONARR(s1) READU,s1unit,sig1 CLOSE,s1unit FREE_LUN,s1unit ; Triangulate over significant parts of images ; triangles0.EU is LONARR(3,115213) ; triangles1.EU is LONARR(3,116598) IF debug THEN PRINT,'Triangulating az-za maps...' OPENR,t0unit,droot+dd+'triangl0.EU',/GET_LUN triangles0=ASSOC(t0unit,LONARR(3,115213)) OPENR,t1unit,droot+dd+'triangl1.EU',/GET_LUN triangles1=ASSOC(t1unit,LONARR(3,116598)) gridspacing=[1.0,1.0] ; 1 deg spacing in az and za boundaries=[0.0,0.0,360.0,72.0] ; minaz,minza,maxaz,maxza ; Generate coefficients for background subtraction coeffs=FLTARR(10,3) ; 10 filters, 3 coeffs per filter IF approx THEN BEGIN coeffs(*,0)=[0.,0.,1.,0.,1.,0.,0.,0.,1.,0.] coeffs(*,1)=[0.,0.,0.,1.,0.,0.,0.,0.,0.,0.] coeffs(*,2)=[1.,0.,0.,0.,0.,1.,0.,1.,0.,0.] same=[0,0,1,1,1,1,0,1,0,0] ENDIF ELSE BEGIN ; Compute Lagrange polynomial coefficients for extrapolating ; backgrounds from measured wavelengths IF debug THEN PRINT,'Generating Lagrange interpolation coefficients...' denom0=(cwl(1)-cwl(9))*(cwl(1)-cwl(6)) denom1=(cwl(9)-cwl(1))*(cwl(9)-cwl(6)) denom2=(cwl(6)-cwl(1))*(cwl(6)-cwl(9)) FOR lw=0,9 DO BEGIN IF (lw NE 1) AND (lw NE 6) AND (lw NE 9) THEN BEGIN coeffs(lw,0)=((cwl(lw)-cwl(9))*(cwl(lw)-cwl(6)))/denom0 coeffs(lw,1)=((cwl(lw)-cwl(1))*(cwl(lw)-cwl(6)))/denom1 coeffs(lw,2)=((cwl(lw)-cwl(1))*(cwl(lw)-cwl(9)))/denom2 ENDIF ; ((lw NE 1) AND (lw NE 6) AND (lw NE 9)) ENDFOR ; (lw) ENDELSE ; (approx) pp=FINDGEN(100)*0.01+0.005 ENDIF ; (bgsubtract) ; read in each cycle of images, decode the headers, rebin to ; (imgsize)x(imgsize) if necessary, scale to 8 bits and display blank=BYTARR(imgsize,imgsize) itime=BYTARR(4) ; load the dark array with representative dark frames for now, ; and supersede them with actual measured ones as they are available darks=INTARR(768,512) ; dummy variables for now rdmeandk,0,1,dk,/nointeractive darks(0,0)=dk rdmeandk,0,10,dk,/nointeractive darks(256,0)=dk rdmeandk,0,60,dk,/nointeractive darks(512,0)=dk rdmeandk,1,1,dk,/nointeractive darks(0,256)=dk rdmeandk,1,10,dk,/nointeractive darks(256,256)=dk rdmeandk,1,60,dk,/nointeractive darks(512,256)=dk dkexp=[1,10,60] procstart=SYSTIME(1) PRINT,procstart-setstart,FORMAT='("Setup took ",F5.1," seconds.")' FOR k=0,N_ELEMENTS(cycle)-1 DO BEGIN ; yield ; This allows IDL to yield control of the processor in ; case another process wants some CPU time. Wait,1 ; Pause execution for 1 sec so IDL can respond to ; keyboard input ; create (or refresh) arrays for temporary storage IF dispim OR (nwsave GT 0) THEN BEGIN ; keep=LONARR(cyclen(k)*256,512) ; ; Solution below implemented 950816 by DPS to avoid ; 'insufficient memory' errors under Windows. IF dos THEN OPENW,keepunit,'keep.tmp',/GET_LUN,/BINARY $ ; IF dos THEN OPENW,keepunit,'keep.tmp',/GET_LUN $ ELSE OPENW,keepunit,'keep.tmp',/GET_LUN keep=ASSOC(keepunit,LONARR(256,256)) ; processed images before background subtraction ; HELP,k,keep keepexst=BYTARR(cyclen(k),2); does image exist? keeptype=keepexst ; image type code, 1-6 keepfilt=keepexst ; filters images taken through keeptime=STRARR(cyclen(k),2); times of images keepexp=keepexst ; exposure lengths keephb=BYTARR(512,2*cyclen(k)) keepaz=FLTARR(cyclen(k),2) ; solar azimuths ENDIF ; (dispim OR (nwsave GT 0)) IF NOT dispim $ THEN PRINT,k,cycle(k),t(cycle(k)),cyclen(k) $ ,FORMAT='("Cycle ",I3," starts at pair "' $ +',I4," (",F6.0," s), and includes "' $ +',I0," pairs")' ip=-1 FOR i=cycle(k),cycle(k)+cyclen(k)-1 DO BEGIN FOR j=0,1 DO BEGIN ip=ip+1 ; update index number for screen position of image p=ip/2 ; position in row ; Is the file path string a valid file name? image_exists=(STRPOS(files(i,j),'___') EQ -1) IF image_exists THEN BEGIN rdkimg,files(i,j),h,im,verbose=debug IF debug THEN PRINT,files(i,j)+' read in' KIh=gethd(h) kihsize=SIZE(KIh) good_header=(kihsize(N_ELEMENTS(kihsize)-2) EQ 8) ENDIF ELSE good_header=0B ; If 'gethd' returned a valid header, carry on. IF good_header THEN BEGIN unitstring=normunits FOR l=0,3 DO itime(l)=KIh.misc.tm.(l) exp=FIX(ROUND(KIh.misc.exp_scale*KIh.exp.exposure/1000.)) exptype=KIh.misc.exposure_type ; If this image is a dark frame, prepare to save it and all others ; in this cycle, otherwise do the appropriate dark subtraction IF exptype EQ 5 THEN BEGIN ; dark frames - save 'em! first_time=(N_ELEMENTS(width) EQ 0) width=512/KIh.exp.sbin ; set up the dark frames length=512/KIh.exp.pbin ; and exp. times arrays IF NOT first_time THEN BEGIN n_darks=N_ELEMENTS(WHERE(dkexp)) old_darks=darks(0:n_darks*width-1,*) old_dkexp=dkexp(0:n_darks-1) ENDIF ; (NOT first_time) ;darks=INTARR(cyclen(k)*width,2*length) ;dkexp=FLTARR(cyclen(k)) IF NOT first_time THEN BEGIN darks(0,0)=old_darks dkexp(0)=old_dkexp ENDIF ; (NOT first_time) IF debug THEN PRINT,'Dark frames being loaded!' dkp=WHERE(dkexp EQ exp,ndkp) ; find exposure IF ndkp EQ 0 $ THEN MESSAGE,'Weird dark exposure!' $ ,/INFORMATIONAL $ ELSE BEGIN ; store the dark frame dkp=dkp(0) darks(dkp*width,j*length)=MEDIAN(im,5) ENDELSE ;IF j EQ 0 THEN dkexp(p)=exp ; and the exposure time ENDIF ELSE BEGIN ; otherwise subtract 'em pd=WHERE(dkexp EQ exp) ; find right exposure time IF pd(0) NE -1 THEN BEGIN pd=pd(0) IF N_ELEMENTS(width) GT 0 THEN BEGIN dwidth=width dlength=length ENDIF ELSE BEGIN dwidth=256 dlength=256 ENDELSE ; (N_ELEMENTS(width) GT 0) dkim=darks(pd*dwidth:(pd+1)*dwidth-1 $ ,j*dlength:(j+1)*dlength-1) ; Ensure the dark frame is the same size as the image. dkcols=(Size(dkim))(1) imcols=(Size(im))(1) CASE 1 OF dkcols GT imcols: dkim=BinF(dkim,dkcols/imcols) dkcols LT imcols: dkim=ReBin(dkim,imcols,imcols,/Sample) $ /(imcols*imcols) ELSE: ENDCASE ;; Following 3 lines added 960822 for debugging purposes ; IF N_Elements(width) EQ 0 THEN BEGIN ; Hlp,dkim,im ; STOP ; ENDIF im=(im-dkim)>0 ENDIF ; (pd(0) NE -1) ENDELSE ; (exptype EQ 5) ; if image isn't 256x256, rebin it IF (KIh.exp.sbin NE 2) OR (KIh.exp.pbin NE 2) THEN BEGIN im=REBIN(im,256,256,/SAMPLE) ; ensure counts are actually summed if a 512x512 image is compressed cfac=512./(256*KIh.exp.sbin) rfac=512./(256*KIh.exp.pbin) im=im*cfac*rfac ENDIF ; ((KIh.exp.sbin NE 2) OR (KIh.exp.pbin NE 2)) ; Apply pseudo-flat-field correction to dark-subtracted image IF (exptype NE 3) AND (exptype NE 5) THEN BEGIN filt=KIh.misc.fw im=ROUND(im*(ffimg(filt*256:(filt+1)*256-1 $ ,j*256:(j+1)*256-1)/100.)) ; Convert to Rayleighs if desired IF rayleighs THEN BEGIN calfac=sens(5*j+filt)/(exp+0.045) IF calfac GT 0 THEN BEGIN im=ROUND(im*calfac) ENDIF ELSE BEGIN unitstring='DN' ENDELSE IF debug THEN PRINT,calfac $ ,FORMAT='("*",F5.2," to get Rayleighs")' IF RperA THEN BEGIN IF (j EQ 0) AND (filt EQ 1) $ THEN im=ROUND(im/rpb(1)) IF (j EQ 1) AND (filt EQ 1) $ THEN im=ROUND(im/rpb(6)) IF (j EQ 1) AND (filt EQ 4) $ THEN im=ROUND(im/rpb(9)) ENDIF ; (RperA) ENDIF ; (rayleighs) ENDIF ; ((exptype NE 3) AND (exptype NE 5)) ; Store processed images for subsequent background subtraction or other use IF dispim OR (nwsave GT 0) THEN BEGIN keepexst(p,j)=1B ; If we're making processed files, despike the images before saving them IF nwsave GT 0 THEN im=cleanmod(im,3,2,3,/quiet) IF debug $ THEN PRINT,files(i,j),pcentile(im,[0.01,0.99]) ; keep(p*256,j*256)=im keep(j*cyclen(k)+p)=im keeptype(p,j)=BYTE(exptype) keepfilt(p,j)=filt keeptime(p,j)=STRING(itime(0:2) $ ,FORMAT='(3(I2.2,:,":"))') keepexp(p,j)=exp sunazel,KIh.misc.dt.year,KIh.misc.dt.month $ ,KIh.misc.dt.day,KIh.misc.tm.hour $ ,KIh.misc.tm.minute,KIh.misc.tm.second $ ,solaz,solel keepaz(p,j)=solaz keephb(0,ip)=h ENDIF ; (dispim OR (nwsave GT 0)) ; Store zenith intensity in appropriate quick-look data file IF ((exptype NE 3) AND (exptype NE 5)) AND $ imsave THEN BEGIN IF j EQ 0 THEN zdata=im(zen0) ELSE zdata=im(zen1) savqldat,j,filt,qlunits,itime,exp,zdata ENDIF ; (((exptype NE 3) AND (exptype NE 5)) AND imsave) ; Save image in quick-look image file for later plotting, if ; appropriate IF KEYWORD_SET(qlimageplot) THEN sto2qlif,unitstring ; Store background images separately for later use IF bgsubtract THEN BEGIN IF (j EQ 0) AND (filt EQ 1) THEN BEGIN bg608=im IF NOT RperA THEN bg608=bg608/rpb(1) ENDIF ; ((j EQ 0) AND (filt EQ 1)) IF (j EQ 1) AND (filt EQ 1) THEN BEGIN bg820=im IF NOT RperA THEN bg820=bg820/rpb(6) ENDIF ; ((j EQ 1) AND (filt EQ 1)) IF (j EQ 1) AND (filt EQ 4) THEN BEGIN bg714=im IF NOT RperA THEN bg714=bg714/rpb(9) ENDIF ; ((j EQ 1) AND (filt EQ 4)) ENDIF ; (bgsubtract) ENDIF ; (valid_header) ENDFOR ; (j) ENDFOR ; (i) FREE_LUN,keepunit ; close ASSOCiated file IF dispim OR (nwsave GT 0) THEN BEGIN OPENR,keepunit,'keep.tmp',/GET_LUN keep=ASSOC(keepunit,LONARR(256,256)) ip=-1 FOR i=cycle(k),cycle(k)+cyclen(k)-1 DO BEGIN ic=i-cycle(k) FOR j=0,1 DO BEGIN ip=ip+1 IF keepexst(ic,j) THEN BEGIN ; identify useful min and max values for byte-scaling ; im=keep(ic*256:(ic+1)*256-1,j*256:(j+1)*256-1) im=keep(j*cyclen(k)+ic) IF dispim THEN BEGIN IF NOT ceiling_set THEN BEGIN IF MAX(im) GT 1 THEN BEGIN IF j EQ 0 $ THEN lmaxi=pcentile(im(sig0),0.995) $ ELSE lmaxi=pcentile(im(sig1),0.995) ENDIF ELSE lmaxi=1 ; (MAX(im) GT 1) ENDIF ELSE lmaxi=topnum(j) ; (NOT ceiling_set) ; if viewing sky image, take dark 'skirt' as min value exptype=keeptype(ic,j) shut=(exptype EQ 3) OR $ ; bias or (exptype EQ 5) ; dark frame IF NOT shut THEN BEGIN IF lmaxi NE 1 THEN BEGIN lmaxbh=pcentile(im(base),0.5)>0 ENDIF ELSE lmaxbh=0 ENDIF ELSE BEGIN ; otherwise take lowest level in dark frame IF j EQ 0 $ THEN lmaxbh=pcentile(im(sig0),0.001) $ ELSE lmaxbh=pcentile(im(sig1),0.001) ENDELSE ; (NOT shut) ; rebin the image to the scale requested by the user IF debug THEN PRINT,lmaxbh,lmaxi,unitstring $ ,FORMAT='("Useful values span "' $ +',I5,"-",I5,1X,A)' IF j EQ 0 THEN mask=mask0 ELSE mask=mask1 IF (((KIh.exp.sbin EQ 2) AND (KIh.exp.pbin EQ 2)) $ AND (fac LE 1.)) $ THEN im=im*mask $ ELSE im=REBIN(im*mask,imgsize,imgsize,/SAMPLE) display(ic*imgsize,j*imgsize)= $ ROTATE(BYTSCL(im,MIN=lmaxbh,MAX=lmaxi $ ,TOP=!D.N_COLORS-1) $ ,5) ENDIF ; (dispim) ; save desired images in associated files fptr=5*j+keepfilt(ic,j) exptype=keeptype(ic,j) IF ((exptype NE 3) AND (exptype NE 5)) $ AND cfsave(fptr) THEN BEGIN hb=keephb(*,ip) ; only save values up to the 99.5th percentile p995=pcentile(im,0.995)>255 ; bbin is the number such that 2**(8+bbin-1) LE p995 LE 2**(8+bbin) bbin=BYTE((CEIL(ALOG((MAX((im>0)0) hb(290)=bbin hb(291)=1B ; save im/(2.**bbin) in a 1-byte space outrec=[hb,REFORM(BYTE(ROUND(im/FLOAT(2^bbin))) $ ,65536,1)] seq=seqnum(fptr) CASE fptr OF 0: save00(seq)=outrec 1: save01(seq)=outrec 2: save02(seq)=outrec 3: save03(seq)=outrec 4: save04(seq)=outrec 5: save10(seq)=outrec 6: save11(seq)=outrec 7: save12(seq)=outrec 8: save13(seq)=outrec 9: save14(seq)=outrec ENDCASE ; (fptr) seqnum(fptr)=seqnum(fptr)+1 ENDIF ; (((exptype NE 3) AND ; (exptype NE 5)) AND cfsave(fptr) ; create and save image annotation strings IF dispim THEN BEGIN camfilt=STRMID(files(i,j),STRLEN(files(i,j))-9,2) cycleparms(0,ip)=camfilt+':'+ext_list(i,j) cycleparms(1,ip)=keeptime(ic,j) cycleparms(2,ip)=STRING(keepexp(ic,j) $ ,FORMAT='(F5.1,1X,"s")') cycleparms(3,ip)=STRING(ROUND(lmaxi),unitstring $ ,FORMAT='("< ",I0.0,A)') cycleparms(4,ip)=STRING(ROUND(lmaxbh),unitstring $ ,FORMAT='("> ",I0.0,A)') ENDIF ; (dispim) ENDIF ELSE BEGIN ; No image found; load blank image and use dummy annotation strings IF dispim THEN BEGIN display(ic*imgsize,j*imgsize)=blank cycleparms(0,ip)='??:-1' cycleparms(1,ip)='??:??:??' cycleparms(2,ip)='??.? s' cycleparms(3,ip)='< 1R' cycleparms(4,ip)='> 0R' ENDIF ; (dispim) ENDELSE ; (keepexst(ic,j)) ENDFOR ; (j) loop through cameras ENDFOR ; (i) loop through images in this cycle ; Blank the rest of the display, if necessary IF dispim THEN BEGIN IF ip LT 2*ninrow-1 THEN BEGIN FOR l=ip+1,2*ninrow-1 DO BEGIN p=l/2 q=l MOD 2 display(p*imgsize,q*imgsize)=blank ENDFOR ; (l) ENDIF ; (ip LT 2*ninrow-1) ENDIF ; (dispim) ; Display the images all at once IF dispim THEN BEGIN cyclestart=cycleparms(1,0) IF STRPOS(cyclestart,'?') NE -1 THEN cyclestart=cycleparms(1,1) PRINT,cyclestart,FORMAT='(/"Ready at ",A,":",$)' IF GET_KBRD(1) EQ 'q' THEN GOTO,DONE mtv,display,res=120. ; Now annotate the images FOR l=0,ip DO BEGIN p=l/2 q=l MOD 2 ; Image camera, filter, and sequence number XYOUTS,scl(res=120.,p*imgsize+1) $ ,scl(res=120.,barwidth+(2-q)*imgsize-X_charheight/fac-1) $ ,cycleparms(0,l) $ ,/DEVICE,ALIGNMENT=0.0,CHARSIZE=csz/fac $ ,COLOR=0.8*!D.N_COLORS ; and image time XYOUTS,scl(res=120.,p*imgsize+1) $ ,scl(res=120.,barwidth+(1-q)*imgsize+X_charheight/fac+3) $ ,cycleparms(1,l) $ ,/DEVICE,ALIGNMENT=0.0,CHARSIZE=csz/fac $ ,COLOR=0.8*!D.N_COLORS ; and exposure time XYOUTS,scl(res=120.,p*imgsize+1) $ ,scl(res=120.,barwidth+(1-q)*imgsize+1) $ ,cycleparms(2,l) $ ,/DEVICE,ALIGNMENT=0.0,CHARSIZE=csz/fac $ ,COLOR=0.8*!D.N_COLORS ; and maximum displayed value XYOUTS,scl(res=120.,(p+1)*imgsize-1) $ ,scl(res=120.,barwidth+(2-q)*imgsize-X_charheight/fac-1) $ ,cycleparms(3,l) $ ,/DEVICE,ALIGNMENT=1.0,CHARSIZE=csz/fac $ ,COLOR=0.8*!D.N_COLORS ; and minimum displayed value XYOUTS,scl(res=120.,(p+1)*imgsize-1) $ ,scl(res=120.,barwidth+(1-q)*imgsize+1) $ ,cycleparms(4,l) $ ,/DEVICE,ALIGNMENT=1.0,CHARSIZE=csz/fac $ ,COLOR=0.8*!D.N_COLORS ; and direction to the Sun solaz=keepaz(p,q) IF solaz EQ 0. THEN sunptr='' ELSE BEGIN IF solaz LT 180. THEN BEGIN sunptr='Sun>' ENDIF ELSE BEGIN sunptr='