PRO getsavhd,filename,nimages,utsec,exptimes,bad_image

ON_IOERROR,go_on

OPENR,hbunit,filename,/GET_LUN
hbfile=ASSOC(hbunit,BYTARR(512))
hbstat=FSTAT(hbunit)
hbbytes=hbstat.SIZE
nimages=hbbytes/66048L
PRINT,'There are '+STRTRIM(nimages,2)+' images'

hbytes=BYTARR(512,nimages)
utsec=LONARR(nimages)
exptimes=FLTARR(nimages)

FOR i=0L,nimages-1 DO BEGIN
    hbytes(0,i)=hbfile(129*i)
    header=gethd(hbytes(*,i))
    utsec(i)=3600L*header.misc.tm.hour $
             + 60L*header.misc.tm.minute $
             +  1L*header.misc.tm.second
    exptimes(i)=header.misc.exp_scale*header.exp.exposure/1000.
    IF (i GT 0) AND (i MOD 20 EQ 0) THEN PRINT,i,FORMAT='(I4,$)'
    IF (i GT 0) AND (i MOD 400 EQ 0) THEN PRINT,FORMAT='()'
    GOTO,endloop

go_on:
    IF KEYWORD_SET(bad_image) $
        THEN bad_image=[bad_image,i] $
        ELSE bad_image=[i]

endloop:
    ENDFOR

PRINT,FORMAT='()'
CLOSE,hbunit
FREE_LUN,hbunit

IF NOT KEYWORD_SET(bad_image) THEN bad_image=[-1]

RETURN
END

;---------------------------------------------------------------------

FUNCTION getlims,nimages,utsec,start=start,finish=finish

lims=INTARR(2)

IF KEYWORD_SET(start) THEN BEGIN
    sstart=SIZE(start)
    IF sstart(0) EQ 0 THEN startut=start $
        ELSE BEGIN
        IF sstart(1) LT 3 THEN start=INTARR(3)+start
        startut=TOTAL([3600.,60.,1.]*start)
        ENDELSE
    after=WHERE(utsec GT startut,nafter)
    IF nafter GT 0 THEN lims(0)=MIN(after) ELSE $
        MESSAGE,'No data after '+STRTRIM(startut,2)+' s UT!'
    ENDIF ELSE lims(0)=0

IF KEYWORD_SET(finish) THEN BEGIN
    sfinish=SIZE(finish)
    IF sfinish(0) EQ 0 THEN finishut=finish $
        ELSE BEGIN
        IF sfinish(1) LT 3 THEN finish=INTARR(3)+finish
        finishut=TOTAL([3600.,60.,1.]*finish)
        ENDELSE
    before=WHERE(utsec LT finishut,nbefore)
    IF nbefore GT 0 THEN lims(1)=MAX(before) ELSE $
        MESSAGE,'No data before '+STRTRIM(finishut,2)+' s UT!'
    ENDIF ELSE lims(1)=nimages-1

RETURN,lims
END

;---------------------------------------------------------------------

PRO getexpt,exptimes,expt,selexp,fexpt=fexpt

; Select the exposure time(s) of the images to be processed.  If
; selexp in non-zero then the value of expt is significant.

mnex=MIN(exptimes,MAX=mxex)
IF KEYWORD_SET(fexpt) THEN BEGIN
    expt=fexpt 
    selexp=1
    ENDIF ELSE BEGIN
    IF mnex NE mxex THEN BEGIN
        REPEAT BEGIN
            PRINT,mnex,mxex $
                 ,FORMAT='("Exposure times range from ",F4.1," to ",F4.1," s")'
            PRINT,'Enter the exposure time of interest (<Return> for all)'
            texp=''
            READ,texp
            selexp=isnumber(texp,expt) 
            j=WHERE(exptimes EQ expt,nj)
            ENDREP UNTIL nj GT 0
        ENDIF ELSE BEGIN
        selexp=0
        expt=FIX(ROUND(mnex))
        ENDELSE
    ENDELSE
RETURN
END

;---------------------------------------------------------------------

FUNCTION gettitle,utsecs,filter,year,month,day,vh,timespan $
                 ,sun=sun,moon=moon,offset=offset $
                 ,van=van,dist=dist
; This routine composes and returns a title based on the parameters
; passed to it.
; FORMAT:
;        'EUREKA (80.0N,86.4W;88.8I'
;        'HH:00-HH:00UT dd MON yyyy'
;        'WWW.W nm Solar azimuth[+DDDd]'
;        'vR hhh km  proj to HHH km'

line0='EUREKA (80.0N,86.4W;88.8!7K!X)'

interval=ROUND(timespan/3600.)        ; hardwired for now, maybe a parameter later
stime=FLOOR(utsecs/timespan)*interval
etime=stime+interval
datestring=ymd2date(year,month,day,format='d$ n$ Y$')
line1=STRING(stime,etime,datestring $
            ,FORMAT='(I2.2,":00-",I2.2,":00UT ",A)')

wavelen=STRMID(filter,0,STRPOS(filter,'nm')+2)
CASE 1 OF
KEYWORD_SET(sun): BEGIN
        az='Solar azimuth'
        IF KEYWORD_SET(offset) THEN $
            IF offset GE 0 THEN $
                az=az+'+'+STRTRIM(offset,2)+'!9%!X' $
                ELSE az=az+'-'+STRTRIM(ABS(offset),2)+'!9%!X'
        END
KEYWORD_SET(moon): BEGIN
        az='Lunar azimuth'
        IF KEYWORD_SET(offset) THEN $
            IF offset GE 0 THEN $
                az=az+'+'+STRTRIM(offset,2)+'!9%!X' $
                ELSE az=az+'-'+STRTRIM(ABS(offset),2)+'!9%!X'
        END
ELSE: BEGIN
        off=0
        IF KEYWORD_SET(offset) THEN off=offset
        az='Azimuth '+STRTRIM(off,2)+'!9%!X'
        END
    ENDCASE
line2=wavelen+' '+az

line3=''
IF KEYWORD_SET(van) THEN $
    line3=line3+'vR '+STRTRIM(ROUND(vh),2)+' km'
IF KEYWORD_SET(dist) THEN $
    line3=line3+'  '+'proj to '+STRTRIM(ROUND(dist),2)+' km'
line3=STRTRIM(line3,2)        ; to lose 2 leading spaces if dist but not van

title=line0+'!C'+line1+'!C'+line2
IF line3 NE '' THEN title=title+'!C'+line3

RETURN,title
END

;---------------------------------------------------------------------

FUNCTION findpeaks,hdist,slice,peak=peak,npeaks=npeaks

; This routine locates peaks in the brightness slice passed from
; the main routine.  The background level, the peak height above
; the background, and the depth of the trough between peaks are
; all settable parameters passed from the main routine.

npeaks=0
IF KEYWORD_SET(peak) THEN BEGIN
    bg=peak(0)
    ht=peak(1)
    depth=peak(2)
    ENDIF ELSE BEGIN
    bg=100
    ht=30
    depth=10
    ENDELSE

IF MAX(slice) GT bg+ht THEN BEGIN
    gslice=slice
    ghdist=hdist
    nvals=N_ELEMENTS(gslice)

; Smooth the profile to a lower spatial resolution
    dh=hdist-SHIFT(hdist,-1)
    bin=10.*CEIL(MAX(dh)/10.)
    hdmin=MIN(hdist,MAX=hdmax)
    dmin=FLOOR(hdmin/bin)*bin
    dmax=FLOOR(hdmax/bin)*bin
    h=HISTOGRAM(hdist,BIN=bin,MIN=dmin,MAX=dmax,REVERSE_INDICES=r)
    nh=N_ELEMENTS(h)
    hdh=FINDGEN(nh)*bin+dmin+bin/2.
    hdm=FLTARR(nh)
    FOR j=0,nh-1 DO BEGIN
        IF r(j) NE r(j+1) THEN BEGIN
            spop=slice(r(r(j):r(j+1)-1))
            hdm(j)=mean(spop)
            ENDIF ELSE hdm(j)=0.
        ENDFOR

; Now find local extrema of binned profile
    blmax=extremes(hdm,+1)
    ; Sort blmax in descending order of peak heights
    blmax=blmax(reverse(SORT(hdm(blmax))))
    blmin=extremes(hdm,-1)
    IF TOTAL(blmin NE -1) EQ 0 THEN mm=MIN(hdm,blmin)
    maxpks=MIN([3,N_ELEMENTS(blmax)])

; Now transfer these to corresponding positions in unbinned slice
    lmax=INTARR(maxpks)
    FOR i=0,maxpks-1 DO BEGIN
        close=WHERE(ABS(ghdist-hdh(blmax(i))) LT bin)
        gpk=MAX(gslice(close),loc)
        lmax(i)=close(0)+loc
        ENDFOR
    lmin=INTARR(N_ELEMENTS(blmin))
    FOR i=0,N_ELEMENTS(lmin)-1 DO BEGIN
        close=WHERE(ABS(ghdist-hdh(blmin(i))) LT bin)
        gtr=MIN(gslice(close),loc)
        lmin(i)=close(0)+loc
        ENDFOR

; Now check whether each local maximum is well-separated from the
; nearest local minimum, and is sufficiently high above the
; neighbouring local minima.  Accept at most 3 maxima per slice.
    keep=BYTARR(maxpks)
    FOR i=0,maxpks-1 DO BEGIN
        IF i EQ 0 THEN keep(i)=(gslice(lmax(i)) GT bg+ht) $
            ELSE BEGIN
            mdist=MIN(ABS(lmax(i)-lmin))
            keep(i)=(gslice(lmax(i)) GT bg+ht) AND (mdist GE 2) 
            ENDELSE
         IF (i GT 0) AND keep(i) THEN BEGIN
            m1=lmax(i-1)
            m2=lmax(i)
            IF m1 GT m2 THEN swap,m1,m2
            locmin=MIN(gslice(m1:m2))
            keep(i)=(locmin LT gslice(lmax(i))-depth)
            ENDIF
        ENDFOR

    IF TOTAL(keep) GE 1 THEN pks=ghdist(lmax(WHERE(keep))) ELSE pks=[-9999]
    ENDIF ELSE pks=[-9999]

npeaks=FIX(TOTAL(pks NE -9999))
RETURN,pks
END

;---------------------------------------------------------------------

PRO aprofile,yr,month,day,cam,filt $
            ,sun=sun,moon=moon,offset=offset,debug=debug,ps=ps $
            ,van=van,vRh=vRh,dist=dist,stack=stack,prompt=prompt $
            ,peak=peak,start=start,finish=finish,fexpt=fexpt
;
; This routine takes the images in a large .RAY file and extracts
; intensities along a meridian from each image.  The profiles can
; then be plotted or transformed to show the progress of emitting
; structures along the meridian.  The meridian can be aligned with
; the Sun (set keyword sun), the Moon (set keyword moon), or offset
; by a fixed angle from either of these reference points.  If neither
; sun nor moon is specified, offset gives the fixed azimuth of the
; meridian.
;
; The keyword 'van' turns on van Rhijn correction, which divides
; the measurements by the shape of a layer of unit vertical brightness
; as seen along a meridian (a van Rhijn profile).  The keyword 'vRh' 
; allows the user to supply an assumed emission height to override 
; the built-in values in the array 'h'.
;
; If a plot of brightness versus distance from overhead rather than
; elevation angle is desired, the keyword 'dist' should be given
; the height (in km) of the layer in which the imaged emission is
; assumed to originate.

; COMMON date,year,month,day

; Initialize
@isitdos
; filters=getfile('/aragorn/user2/steele/idl/filters.doc')
s=''
psplot=KEYWORD_SET(ps)
!ORDER=1
IF !D.WINDOW EQ 0 THEN WDELETE,0
!X.MARGIN=[10,8]
ytop=7
!Y.MARGIN=[4,ytop]

; ensure a 4-digit year for display; make reasonable assumptions
year=yr
IF year LT 1900 THEN BEGIN
    IF year LT 50 THEN year=year+2000
    IF year GT 50 THEN year=year+1900
ENDIF

IF KEYWORD_SET(van) THEN $
    h=[87.,100.,300.,150.,300.,87.,100.,150.,97.,100.]

; Find out which pixels are viewing the sky
getsig,year,month,day,sig0,sig1
CASE cam OF
 0: sig=sig0
 1: sig=sig1
    ENDCASE

; Get the path to the desired images
template=FiLocate(year,month,day,cam,filt)

; Open the file as a series of small (header-size) units and get
; the header information
getsavhd,template,nimages,utsec,exptimes,bad_image

; Identify first and last image to process
lims=getlims(nimages,utsec,start=start,finish=finish)

; Set up azimuth and zenith angle images for the specified camera
; and site
PCsite=PoCasite(year,month,day)
col=BINDGEN(256)#REPLICATE(1B,1,256)
row=TRANSPOSE(col)
cr2az,year,month,day,cam,filt,col,row,az,za,site=PCsite,/degrees
IF KEYWORD_SET(offset) THEN offset=FIX(ROUND(offset))

; identify image sequences and prompt user for sequence(s) to display
dt=ABS(utsec-SHIFT(utsec,1))        ; rectify big negative number in front
seqstart=WHERE(dt GT 600)        ; semi-arbitrary large number (10 min)

;; Set plotting device
;IF NOT KEYWORD_SET(ps) THEN BEGIN
;    resp=''
;    PRINT,'Make a PostScript plot file? (Y/[N])'
;    READ,resp
;    psplot=(STRUPCASE(resp) EQ 'Y')
;    ENDIF

IF psplot THEN BEGIN
    psname=STRING(year MOD 100,month,day,cam,filt $
                 ,FORMAT='(3I2.2,2I1,".sps")')
    psopen,/long,psname
    ENDIF

IF seqstart(0) GE 0 THEN BEGIN
    PRINT,'Sequences begin with the following image numbers(s):'
    PRINT,seqstart
    ENDIF

; Select exposure time(s) of images to be processed
getexpt,exptimes,expt,selexp,fexpt=fexpt

; Open the single big file and associate it to variable.
OPENR,imunit,template,/GET_LUN
files=ASSOC(imunit,BYTARR(66048))

; If a stack plot is desired, set up the big window and the axes
IF KEYWORD_SET(stack) THEN BEGIN
    IF !D.NAME NE 'PS' THEN BEGIN
        ss=Screen_Size()
        WINDOW,0,XSIZE=0.64*ss(1),YSIZE=ss(1)
        ENDIF
    !P.POSITION=[0.0,0.0,0.90,0.98]
    !X.OMARGIN=[8,8]
    EMPTY
    !P.MULTI=[0,2,1,0,0]
    ENDIF

; Initialize quantities needed below
ip=-1                                ; counter for images actually processed
olduts=86399.                        ; dummy - no longer needed
fullscale=8000.                        ; limit of intensity scale for stack plot
fext=(fullscale+250.)/fullscale        ; we allow 250 R at top of each graph
timespan=14400.                        ; duration (s) of each stack plot
utnewpage=-2*timespan-1.        ; dummy time of last axis start, superseded below
ftimespan=FIX(ROUND(timespan))        ; sometimes need an integer
nticks=ROUND(timespan/900.)        ; # of tick intervals needed for 15-min resolution
tvals=FINDGEN(nticks+1)                ; times (s) for use in plotting time axis

; read in each image, decode the header, rescale from 8 bits and process
FOR i=lims(0),lims(1)-1 DO BEGIN
    whbad=WHERE(bad_image EQ i,nwhbad)
    IF nwhbad EQ 0 THEN BEGIN
        file=files(i)
        im=BYTE(file,512,256,256)
    ; If file was renormalized to fit 8 bits, reverse this.
        IF file(290) GT 0 THEN BEGIN
            IF file(291) EQ 1 THEN im=LONG(im)*(2^file(290))
            IF file(291) EQ 0 THEN im=LONG(im)*file(290)
            ENDIF
        KIh=gethd(BYTE(file,0,512))

; extract exposure time if necessary
        IF selexp THEN BEGIN
            thisexp=KIh.misc.exp_scale*KIh.exp.exposure/1000.
            go_on=(ABS(thisexp-expt) LT 1.0)
            ENDIF ELSE go_on=1

        IF go_on THEN BEGIN
; extract start time of image from header
            itime=BYTARR(4)
            FOR itptr=0,3 DO itime(itptr)=KIh.misc.tm.(itptr)
            utsecs=TOTAL(itime*[3600,60,1,.01])

; Calculate desired azimuth for this image
            CASE 1 OF
 KEYWORD_SET(sun): 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 $
                          ,meraz
KEYWORD_SET(moon): moonazel,KIh.misc.dt.year,KIh.misc.dt.month,KIh.misc.dt.day $
                           ,KIh.misc.tm.hour,KIh.misc.tm.minute,KIh.misc.tm.second $
                           ,meraz
             ELSE: meraz=0.
                 ENDCASE
            IF KEYWORD_SET(offset) THEN meraz=meraz+offset
            meraz=meraz MOD 360

            ip=ip+1                ; index number for screen position of image

            merslice,im*sig,az,za,meraz,elev,slice,sd $
                    ,swidth=1.0,debug=debug

;            IF NOT KEYWORD_SET(stack) THEN BEGIN
                mertitle=STRTRIM(ROUND(meraz),2)+' !9%!X azimuth'
                IF KEYWORD_SET(sun) THEN BEGIN
                    mertitle='Solar Azimuth ('+STRTRIM(ROUND(meraz),2)+' !9%!X)'
                    IF KEYWORD_SET(offset) THEN $
                        mertitle=STRTRIM(offset,2)+' !9%!X offset from '+mertitle
                    ENDIF ELSE IF KEYWORD_SET(moon) THEN BEGIN
                        mertitle='Lunar Azimuth ('+STRTRIM(ROUND(meraz),2)+' !9%!X)'
                        IF KEYWORD_SET(offset) THEN $
                            mertitle=STRTRIM(offset,2)+' !9%!X offset from '+mertitle
                        ENDIF
                mertitle=STRING(itime(0),itime(1),itime(2),FORMAT='(3(I2.2,:,":"))') $
                         +'  '+mertitle
;                ENDIF ELSE BEGIN

;                ENDELSE

; Is van Rhijn correction desired?
            vh=0.
            IF KEYWORD_SET(van) THEN BEGIN
                IF KEYWORD_SET(vRh) THEN vh=vRh ELSE vh=h(5*cam+filt)
                v=vanrhijn(ABS(elev-90.),vh)
                slice=slice/v
                ENDIF

; Exclude points lying outside active image area of CCD
            wmer=WHERE(sd)
            elev=elev(wmer)
            slice=slice(wmer)
            sd=sd(wmer)

; Try median smoothing to lose single-pixel spikes
            slice=MEDIAN(slice,3)

; Brightness versus distance desired
            IF KEYWORD_SET(dist) THEN BEGIN
                eflag=1
                eheight=dist
                hdist=za2dist(elev-90.,eheight,debug=debug,rad=re)

; Smoothing works well.  Now try identifying peaks and troughs
; where background level, peak amplitude, and interpeak trough
; depth are all settable parameters: peak=[back,height,depth]
                IF KEYWORD_SET(peak) THEN $
                    pks=findpeaks(hdist,slice,npeaks=npeaks,peak=peak) $
                    ELSE npeaks=0
                ENDIF

            IF NOT KEYWORD_SET(ymax) THEN ymax=nicenumber(2*MAX(slice),/floor)
            IF NOT psplot THEN WSET,0
            arg=KEYWORD_SET(prompt)
            IF STRUPCASE(GET_KBRD(arg)) EQ 'Q' THEN GOTO,DONE

; Plot versus desired independent variable 
            CASE 1 OF
KEYWORD_SET(dist): BEGIN
                    IF ip EQ 0 THEN maxdist=50*CEIL((2.75*dist)/50)
                    IF KEYWORD_SET(stack) THEN BEGIN
                        toff=(utsecs MOD ftimespan)*(fullscale/timespan)
                        nvals=N_ELEMENTS(hdist)
                        left=MAX([maxdist,hdist+1])
                        right=MIN([-maxdist,hdist-1])
                        hdist=[left,hdist(0),hdist,hdist(nvals-1),right]
                        slice=[0.,0.,slice,0.,0.]
                        PRINT,itime(0),itime(1),itime(2) $
                             ,FORMAT='(I2.2,":",I2.2,":",I2.2," ",$)'
                        IF ((ip+1) MOD 11) EQ 0 THEN PRINT,FORMAT='()'
                        CASE 1 OF

; Time to start a new page of plotting
;;;;(utsecs MOD (2*ftimespan)) LT (olduts MOD (2*ftimespan)): BEGIN
(FLOOR(utsecs)/(2*ftimespan) GT FLOOR(utnewpage)/(2*ftimespan)) OR $
                 (ip EQ 0): BEGIN
                                IF ip EQ 0 THEN multiplot $
                                    ELSE BEGIN
                                    ERASE
                                    multiplot,/reset
                                    multiplot
                                    ENDELSE
                                ptitle=gettitle(utsecs,filters(5*cam+filt) $
                                               ,year,month,day,vh,timespan $
                                               ,sun=sun,moon=moon,offset=offset $
                                               ,van=van,dist=dist)
                                PLOT,hdist,slice+toff $
                                    ,XTITLE='Dist. toward specified azimuth (km)' $
                                    ,YTITLE='Intensity (R)' $
                                    ,TITLE=ptitle $
                                    ,XRANGE=[maxdist,-maxdist],/XSTYLE $
                                    ,YRANGE=[0,fullscale+250],/YSTYLE $
                                    ,XTICKLEN=0.01,YTICKLEN=0.04
                                utnewpage=utsecs
                                utnewaxis=utsecs
                                IF npeaks GT 0 THEN $
                                    FOR k=0,npeaks-1 DO BEGIN
                                        tmp=interpol(slice,hdist,pks(k))
                                        PLOTS,[pks(k),pks(k)] $
                                             ,[toff,toff+tmp] $
                                             ,LINESTYLE=0
                                        PLOTS,[pks(k)] $
                                             ,[toff] $
                                             ,PSYM=1
                                        ENDFOR
                                END

; Time to start the second half of the present page of plotting
FLOOR(utsecs)/ftimespan GT FLOOR(utnewaxis)/ftimespan: BEGIN
;(utsecs MOD ftimespan) LT (olduts MOD ftimespan): BEGIN
                                multiplot
                                ptitle=gettitle(utsecs,filters(5*cam+filt) $
                                               ,year,month,day,vh,timespan $
                                               ,sun=sun,moon=moon,offset=offset $
                                               ,van=van,dist=dist)
                                PLOT,hdist,slice+toff $
                                    ,XTITLE='Dist. toward specified azimuth (km)' $
                                    ,TITLE=ptitle $
                                    ,XRANGE=[maxdist,-maxdist],/XSTYLE $
                                    ,YRANGE=[0,fullscale+250],/YSTYLE $
                                    ,XTICKLEN=0.01,YTICKLEN=0.04
                                utnewaxis=utsecs
                                IF npeaks GT 0 THEN $
                                    FOR k=0,npeaks-1 DO BEGIN
                                        tmp=interpol(slice,hdist,pks(k))
                                        PLOTS,[pks(k),pks(k)] $
                                             ,[toff,toff+tmp] $
                                             ,LINESTYLE=0
                                        PLOTS,[pks(k)] $
                                             ,[toff] $
                                             ,PSYM=1
                                        ENDFOR
                                AXIS,-maxdist,YAXIS=1 $
                                    ,YTITLE='Time within plotted interval (HH:MM)' $
                                    ,YTICKLEN=-!P.TICKLEN $
                                    ,/DATA,YTICKS=nticks $
                                    ,YTICKV=tvals*(fullscale/nticks) $
                                    ,YTICKNAME=time_label(tvals*900)
                                tit_ht=CONVERT_COORD(maxdist,fullscale+250.,0.  $
                                                   ,/DATA,/TO_DEVICE)
                                PLOTS,[tit_ht(0),tit_ht(0)] $
                                     ,[tit_ht(1),tit_ht(1)+(ytop-0.5)*!D.Y_CH_SIZE] $
                                     ,/DEVICE
                                END

; Continue plotting on the present axes
                      ELSE: BEGIN
                                OPLOT,hdist,slice+toff
                                IF npeaks GT 0 THEN $
                                    FOR k=0,npeaks-1 DO BEGIN
                                        tmp=interpol(slice,hdist,pks(k))
                                        PLOTS,[pks(k),pks(k)] $
                                             ,[toff,toff+tmp] $
                                             ,LINESTYLE=0
                                        PLOTS,[pks(k)] $
                                             ,[toff] $
                                             ,PSYM=1
                                        ENDFOR
                                END
                            ENDCASE
                        ENDIF ELSE BEGIN
                        PLOT,hdist,slice,TITLE=mertitle $
                            ,XTITLE='Dist. toward specified azimuth (km)' $
                            ,YTITLE='Intensity (R)' $
                    ;            ,PSYM=3 $
                            ,XRANGE=[maxdist,-maxdist],/XSTYLE $
                            ,YRANGE=[0,ymax],/YSTYLE
                        OPLOT,hdist,sd,COLOR=!D.N_COLORS/2
                        ENDELSE
;                    XYOUTS,0.03*!X.CRANGE(0)+0.97*!X.CRANGE(1) $
;                          ,0.05*!Y.CRANGE(0)+0.95*!Y.CRANGE(1) $
;                          ,'Emission height '+STRTRIM(dist,2)+' km' $
;                          ,ALIGNMENT=1
                    END
          ELSE: BEGIN
                    PLOT,elev,slice,TITLE=mertitle $
                        ,XTITLE='Elevation angle (!9%!X)' $
                        ,YTITLE='Intensity (R)' $
                ;        ,PSYM=3 $
                        ,XTICKS=6,XTICKV=FINDGEN(7)*30 $
                        ,XRANGE=[0,180],/XSTYLE $
                        ,YRANGE=[0,ymax],/YSTYLE
                    OPLOT,elev,sd,COLOR=!D.N_COLORS/2
                    END
                ENDCASE

;            IF KEYWORD_SET(van) THEN $
;                XYOUTS,0.97*!X.CRANGE(0)+0.03*!X.CRANGE(1) $
;                      ,0.05*!Y.CRANGE(0)+0.95*!Y.CRANGE(1) $
;                      ,'van Rhijn corrected'
            EMPTY

            ENDIF
        ENDIF
    ENDFOR

DONE:
; Ensure a time axis is shown on the last page even if only one
; panel is plotted
AXIS,-maxdist,YAXIS=1 $
    ,YTITLE='Time within plotted interval (HH:MM)' $
    ,YTICKLEN=-!P.TICKLEN $
    ,/DATA,YTICKS=nticks $
    ,YTICKV=tvals*(fullscale/nticks) $
    ,YTICKNAME=time_label(tvals*900)

IF psplot THEN BEGIN
    PRINT,FORMAT='()'
    psclose
;    SPAWN,'~steele/scripts/rmps3 '+psname+' b'
    ENDIF ELSE BEGIN
    PRINT,FORMAT='()'
    PRINT,STRING(7B)+'Press any key to erase the display and end.'
    j=GET_KBRD(1)
    WDELETE,2
    WDELETE,0
    ENDELSE

END