PRO dmspct,tcol,bw=bw

; If keyword bw is set, use a grey scale colour table.

IF KEYWORD_SET(bw) THEN BEGIN
    r=INDGEN(15)/14.
    g=r
    b=r
    tcol=1B

; Establish the number of colors and the relative RGB intensities.

    ENDIF ELSE BEGIN
    r=[0.2,0.5,0.2,0.5,0.5,0.3,0.0,0.0,0.6,1.0,1.0,1.0,1.0,1.0,1.0]
    g=[0.0,0.0,0.0,0.3,0.5,0.6,0.5,0.8,0.9,1.0,0.8,0.5,0.0,0.5,0.8]
    b=[0.2,0.5,0.6,0.9,0.9,0.9,0.0,0.0,0.0,0.0,0.2,0.2,0.0,0.8,0.8]
    tcol=0.35*!D.N_COLORS
    ENDELSE

; Expand to !D.N_COLORS elements and normalize to 256 steps.

samp=BYTE(FINDGEN(!D.N_COLORS)*(15./!D.N_COLORS))
rexp=BYTSCL(r(samp))
gexp=BYTSCL(g(samp))
bexp=BYTSCL(b(samp))

; Load the color table.

TVLCT,rexp,gexp,bexp

RETURN
END

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

FUNCTION power10,axis,index,value

power=ALOG10(value)
nice=(ABS(power-ROUND(power)) LT 0.01)
IF nice $
    THEN powst=STRTRIM(ROUND(power),2) $
    ELSE powst=STRTRIM(STRING(power,FORMAT='(F12.2)'),2)
RETURN,'10!e'+powst+'!n'
END

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

FUNCTION nullst,axis,index,value

RETURN,' '
END

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

;FUNCTION GEOMF(I,J)
FUNCTION geomf,j
;
;       NOTE: The single parameter to the IDL function gives the
;	satellite number (6..12).
;
;                          ADDED TO CODE 25 AUG 1992
;
;     I = CHANNELS 1 TO 40
;     J = ('F6', 'F7', 'F8', 'F9', 'F10', 'F11' AND 'F12')
;
;     GEOM RETURNS THE GEOMETRIC FACTOR WITH UNITS (CM**2*STER*KEV)
;
;	Converted from FORTRAN to IDL by DPS on 950519.
;
;     SAVE
COMMON geomfact,geof6,geof7,geof8,geof9,geof10,geof11,geof12
;
;     CHARACTER*3 J
;     DIMENSION GEOF6(40),GEOF7(40),GEOF8(40),GEOF9(40),GEOF10(40)
;    +         ,GEOF11(40),GEOF12(40)
;
geof6=[75.0,49.0,41.0,33.0,27.0,21.0,16.0,13.0,9.60, $
       7.60,1.40,1.00,0.80,0.51,0.34,0.20,0.10,0.05, $
       0.02,.0063,180.,120.,81.0,55.0,38.0,26.0,18.0, $
       12.0,8.10,5.60,230.,150.,140.,79.0,55.0,37.0, $
       25.0,18.0,12.0,8.00]

geof7=[58.0,49.0,41.0,33.0,27.0,21.0,16.0,13.0,9.60, $
       7.60,3.30,2.50,2.00,1.30,0.91,0.58,0.34,0.20, $
       0.11,.057,240.,160.,110.,75.0,52.0,35.0,24.0, $
       17.0,11.0,7.60,210.,140.,98.0,66.0,46.0,31.0, $
       21.0,15.0,10.0,7.00]

geof8=[32.6,27.5,23.0,18.5,15.2,11.8,8.99,7.30,5.39, $
       4.27,3.20,2.43,1.94,1.26,.884,.563,.330,.194, $
       .107,.0553,115.,76.9,52.9,36.1,25.0,16.8,11.5, $
       8.17,5.29,3.65,111.,74.1,51.9,34.9,24.3,16.4, $
       11.1,7.94,5.29,3.70]

geof9=[20.1,17.0,14.2,11.5,9.37,7.29,5.55,4.51,3.33, $
       2.64,3.93,2.98,2.38,1.55,1.08,0.69,.405,.238, $
       .131,.0678,114.,76.2,52.4,35.7,24.8,16.7,11.4, $
       8.10,5.24,3.62,130.,86.4,60.5,40.7,28.4,19.1, $
       13.0,9.26,6.17,4.32]

geof10=[46.2,38.9,30.2,24.3,19.4,15.8,12.6,10.2,8.24, $
        6.69,6.37,4.32,2.92,1.99,1.37,.930,.637,.438, $
        .302,.206,56.7,38.7,26.0,17.5,11.9,8.18,5.54, $
        3.78,2.57,1.77,50.1,33.9,23.0,15.7,10.8,7.31, $
        5.01,3.45,2.38,1.62]

geof11=[63.1,54.4,42.8,34.9,27.8,22.6,18.2,14.9,12.1, $
        9.50,10.1,6.85,4.71,3.10,2.09,1.43,.963,.663, $
        .449,.300,54.9,38.1,26.0,17.7,12.1,8.24,5.64, $
        3.87,2.66,1.77,50.1,35.0,24.0,16.7,11.4,7.83, $
        5.27,3.71,2.56,1.77]

geof12=[41.4,34.7,29.0,24.2,20.3,16.9,14.2,11.8,9.90, $
        8.28,7.30,4.91,3.30,2.22,1.49,1.00,.673,.453, $
        .304,.204,70.6,48.2,32.9,22.5,15.3,10.5,7.15, $
        4.88,3.33,2.28,40.2,27.5,18.8,12.9,8.82,6.04, $
        4.13,2.83,1.94,1.32]

      geofac=REPLICATE(-1.,40)
;     IF(J.EQ.' F6') THEN
;        GEOMF=GEOF6(I)
;     ELSE IF(J.EQ.' F7') THEN
;        GEOMF=GEOF7(I)
;     ELSE IF(J.EQ.' F8') THEN
;        GEOMF=GEOF8(I)
;     ELSE IF(J.EQ.' F9') THEN
;        GEOMF=GEOF9(I)
;     ELSE IF(J.EQ.'F10') THEN
;        GEOMF=GEOF10(I)
;     ELSE IF(J.EQ.'F11') THEN
;        GEOMF=GEOF11(I)
;     ELSE IF(J.EQ.'F12') THEN
;        GEOMF=GEOF12(I)
;     ENDIF
CASE j OF
 6: geofac=geof6
 7: geofac=geof7
 8: geofac=geof8
 9: geofac=geof9
10: geofac=geof10
11: geofac=geof11
12: geofac=geof12
ELSE: MESSAGE,'Bad satellite: '+STRTRIM(j,2)
    ENDCASE
;
IF TOTAL(geofac) GT 0. THEN geofac=geofac/100000.
;
RETURN,geofac
END

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

;      FUNCTION DIFF(K,J)
FUNCTION diff,j
;
;       NOTE: The single parameter to the IDL function is the
;	satellite number 6..12.
;
;C                         ADDED TO CODE:  25 AUG 1992
;      SAVE
;      CHARACTER*3 J
;      DIMENSION DIFF6(40),DIFF7(40),DIFF8(40),DIFF9(40)
;      DIMENSION DIFF10(40),DIFF11(40),DIFF12(40)
;C
;      DATA NCALL/0/
COMMON diffstuff,ncall,diff6,diff7,diff8,diff9,diff10,diff11,diff12
;C
;      IF(NCALL.EQ.0) THEN
IF N_ELEMENTS(diff6) EQ 0 THEN BEGIN
    ncall=1B
;	NCALL=1
;	DO 10 I=1,40
;	  DIFF6(I)=1./(GEOMF(I,' F6')*.098)
    diff6=1./(geomf(6)*0.098)
;	  DIFF7(I)=1./(GEOMF(I,' F7')*.098)
    diff7=1./(geomf(7)*0.098)
;	  DIFF8(I)=1./(GEOMF(I,' F8')*.098)
    diff8=1./(geomf(8)*0.098)
;	  DIFF9(I)=1./(GEOMF(I,' F9')*.098)
    diff9=1./(geomf(9)*0.098)
;	  DIFF10(I)=1./(GEOMF(I,'F10')*.098)
    diff10=1./(geomf(10)*0.098)
;	  DIFF11(I)=1./(GEOMF(I,'F11')*.098)
    diff11=1./(geomf(11)*0.098)
;	  DIFF12(I)=1./(GEOMF(I,'F12')*.098)
    diff12=1./(geomf(12)*0.098)
;  10    CONTINUE
;      ENDIF
    ENDIF
;C
;      DIFF=-1
diffp=REPLICATE(-1.,40)
;      IF(K.LT.1.OR.K.GT.40) RETURN
;C
;      IF(J.EQ.' F6') THEN
;	DIFF=DIFF6(K)
;      ELSE IF(J.EQ.' F7') THEN
;	DIFF=DIFF7(K)
;      ELSE IF(J.EQ.' F8') THEN
;	DIFF=DIFF8(K)
;      ELSE IF(J.EQ.' F9') THEN
;	DIFF=DIFF9(K)
;      ELSE IF(J.EQ.'F10') THEN
;	DIFF=DIFF10(K)
;      ELSE IF(J.EQ.'F11') THEN
;	DIFF=DIFF11(K)
;      ELSE IF(J.EQ.'F12') THEN
;	DIFF=DIFF12(K)
;      ENDIF
CASE j OF
 6: diffp=diff6
 7: diffp=diff7
 8: diffp=diff8
 9: diffp=diff9
10: diffp=diff10
11: diffp=diff11
12: diffp=diff12
ELSE:
    ENDCASE
;C
;      RETURN
;      END
RETURN,diffp
END

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

;      FUNCTION ENG(K,I)
FUNCTION eng,k
;
;       NOTE: The single parameter to the IDL function is the
;	satellite number 6..12.
;
;      SAVE
;      CHARACTER*3 I
;      DIMENSION EN6(20),EN7(20),EN8(20),EN9(20),ENG10(20),ENG11(40)
COMMON engstuff,en6,en7,en8,en9,eng10,eng11
;C
;      DATA EN6/30180.,20620.,14040.,9580.,6500.,4420.,3050.,2060.,1410.,
;     +  984.,992.,679.,462.,317.,213.,145.,100.,68.,46.,32./
en6=[30180.,20620.,14040.,9580.,6500.,4420.,3050.,2060.,1410., $
     984.,992.,679.,462.,317.,213.,145.,100.,68.,46.,32.]
;C
;      DATA EN7/30340.,20710.,14070.,9570.,6570.,4470.,3050.,2100.,1430.,
;     +         985.,995.,679.,462.,314.,215.,147.,100.,69.,47.,32./
en7=[30340.,20710.,14070.,9570.,6570.,4470.,3050.,2100.,1430., $
     985.,995.,679.,462.,314.,215.,147.,100.,69.,47.,32.]
;C
;      DATA EN8/31300.,21100.,14300.,9720.,6610.,4500.,3050.,2070.,1400.,
;     +         950.,950.,640.,440.,310.,210.,144.,98.,68.,45.,31./
en8=[31300.,21100.,14300.,9720.,6610.,4500.,3050.,2070.,1400., $
     950.,950.,640.,440.,310.,210.,144.,98.,68.,45.,31.]
;C
;      DATA EN9/31300.,21100.,14300.,9720.,6610.,4500.,3050.,2070.,1400.,
;     +         950.,950.,640.,440.,310.,210.,144.,98.,68.,45.,31./
en9=[31300.,21100.,14300.,9720.,6610.,4500.,3050.,2070.,1400., $
     950.,950.,640.,440.,310.,210.,144.,98.,68.,45.,31.]
;C
;      DATA ENG10/30500.,20800.,14000.,9400.,6400.,4400.,2980.,2030.,
;     +           1380.,950.,960.,650.,440.,300.,206.,140.,96.,66.,46.,
;     +           31./
eng10=[30500.,20800.,14000.,9400.,6400.,4400.,2980.,2030., $
       1380.,950.,960.,650.,440.,300.,206.,140.,96.,66.,46., $
       31.]
;C
;      DATA ENG11/29500.,20500.,14000.,9500.,6500.,4430.,3030.,2080.,
;     +           1430.,950.,960.,670.,460.,320.,218.,150.,101.,71.,49.,
;     +           34.,29500.,20500.,14000.,9500.,6500.,4430.,3030.,2080.,
;     +           1430.,950.,940.,640.,440.,290.,195.,134.,90.,62.,42.,
;     +           28./
eng11=[29500.,20500.,14000.,9500.,6500.,4430.,3030.,2080., $
       1430.,950.,960.,670.,460.,320.,218.,150.,101.,71.,49., $
       34.,29500.,20500.,14000.,9500.,6500.,4430.,3030.,2080., $
       1430.,950.,940.,640.,440.,290.,195.,134.,90.,62.,42., $
       28.]
;C
;      IF(I.EQ.' F6')THEN
;	ENG=EN6(MOD(K-1,20)+1)
;      ELSE IF(I.EQ.' F8') THEN
;	ENG=EN8(MOD(K-1,20)+1)
;      ELSE IF(I.EQ.' F9') THEN
;	ENG=EN9(MOD(K-1,20)+1)
;      ELSE IF(I.EQ.' F7') THEN
;	ENG=EN7(MOD(K-1,20)+1)
;      ELSE IF(I.EQ.'F10') THEN
;	ENG=ENG10(MOD(K-1,20)+1)
;      ELSE
;	ENG=ENG11(K)
;      ENDIF
CASE k OF
 6: engp=[en6,en6]
 7: engp=[en7,en7]
 8: engp=[en8,en8]
 9: engp=[en9,en9]
10: engp=[eng10,eng10]
11: engp=eng11
12: engp=eng11		; assume F12 is similar to F11
ELSE:
    ENDCASE
;C
;      ENG=ENG/1000.
engp=engp/1000.
;C
;      RETURN
;      END
RETURN,engp
END

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

PRO getstep,sutsec,eutsec,utstep,psutsec,peutsec,debug=debug

utopt=[1.,5.,10.,15.,20.,30.,60.,90.,120.,150.,180.,240.,300.]
nopt=N_ELEMENTS(utopt)
div=10.			; this divisor may need adjusting

utspan=eutsec-sutsec
IF KEYWORD_SET(debug) THEN HELP,utspan
utseg=utspan/div
IF KEYWORD_SET(debug) THEN HELP,utseg

poss=WHERE(utopt GE utseg,nposs)
IF KEYWORD_SET(debug) THEN PRINT,'poss:',FIX(poss)
best=poss(0)
IF KEYWORD_SET(debug) THEN HELP,best
utstep=utopt(best)	
IF KEYWORD_SET(debug) THEN HELP,utstep

psutsec=utstep*FLOOR(sutsec/utstep)
IF KEYWORD_SET(debug) THEN HELP,psutsec
peutsec=utstep*CEIL(eutsec/utstep)
IF KEYWORD_SET(debug) THEN HELP,peutsec

WHILE ((peutsec-psutsec)/utstep GT div) AND (best LT (nopt-1)) DO BEGIN
    best=(best+1)<(nopt-1)
    IF KEYWORD_SET(debug) THEN HELP,best
    utstep=utopt(best)
    IF KEYWORD_SET(debug) THEN HELP,utstep
    psutsec=utstep*FLOOR(sutsec/utstep)
    IF KEYWORD_SET(debug) THEN HELP,psutsec
    peutsec=utstep*CEIL(eutsec/utstep)
    IF KEYWORD_SET(debug) THEN HELP,peutsec
    ENDWHILE

RETURN
END

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

PRO rdj4,name,stime,etime $
	,debug=debug,out=out,bw=bw,fpalt=fpalt,spectra=spectra

; This is starting out as a halting attempt to make sense of raw
; DMSP particle data acquired from AFGL.  The VAX/VMS RMS files
; have fixed 5280-byte records, different from the standard 7560-byte
; binary data from the NOAA archive.  A FORTRAN program in Ed
; Weber's directory gives a little hint about how to make sense of
; the data.  This code tries to take it from there, with major help
; from Bob Sheehan at PL.

COMMON geomfact,geof6,geof7,geof8,geof9,geof10,geof11,geof12
COMMON diffstuff,ncall,diff6,diff7,diff8,diff9,diff10,diff11,diff12
COMMON engstuff,en6,en7,en8,en9,eng10,eng11

: Set up OS-dependent parameters
@isitdos

IF N_PARAMS() EQ 0 THEN BEGIN
    MESSAGE,'Usage: rdj4,name,start_hhmm,finish_hhmm[,/debug,/out,/bw,fpalt=fpalt]' $
	   ,/INFORMATIONAL
    RETURN
    ENDIF

name=j4root+dd+name
f=rfindfile(name,COUNT=nf)
IF nf EQ 0 THEN BEGIN
    MESSAGE,'File not found: '+name,/INFORMATIONAL
    RETURN
    ENDIF ELSE BEGIN

    ; Figure out which satellite the data are from

    uname=STRUPCASE(name)
    IF dos THEN BEGIN
        satnum=STRPOS(uname,1,1)
        IF NOT isnumber(satnum) THEN BEGIN	; assume hex sat. num.
            offset=BYTE(satnum)-BYTE('A')
            satnum=10+offset
            ENDIF
        ENDIF ELSE BEGIN
        satpos=STRPOS(uname,'F')
        IF satpos GE 0 THEN BEGIN
	    satnum=FIX(STRMID(uname,satpos+1,1)) 
	    IF satnum LT 6 $	; F10 or later
                THEN satnum=10*satnum+FIX(STRMID(uname,satpos+2,1))
	    ENDIF ELSE BEGIN
	    READ,'Enter the satellite number, please: (e.g., 8 for F8) ' $
                ,satnum
            ENDELSE
        ENDELSE

    ; Get the file size and # of records

    OPENR,lun,f(0),/GET_LUN;,RECORD_LENGTH=5280
    stat=FSTAT(lun)
    flen=stat.SIZE
    nrec=flen/5280
    IF KEYWORD_SET(debug) THEN PRINT,nrec

    ; Load the headers first

    headers=INTARR(15,nrec)
    hb=ASSOC(lun,INTARR(2640))
    FOR i=0,nrec-1 DO BEGIN
	temphead=hb(0:14,i)
	IF !VERSION.OS EQ 'RISC/os' THEN BYTEORDER,temphead
	headers(0,i)=temphead
	ENDFOR
    FREE_LUN,lun

    ; Decode the headers

    headers(4,*)=headers(4,*)+50
    headers(5,*)=headers(5,*)-900
    headers(8,*)=headers(8,*)-900
    headers(10,*)=headers(10,*)-900
    ; The following lines recover southern hemisphere latitudes
    ws=WHERE(headers(5,*) GE '800'XL,nws)
    IF nws GT 0 THEN headers(5,ws)=(headers(5,ws) OR 'FFFFF000'XL)+1
    ws=WHERE(headers(8,*) GE '800'XL,nws)
    IF nws GT 0 THEN headers(8,ws)=(headers(8,ws) OR 'FFFFF000'XL)+1
    ws=WHERE(headers(10,*) GE '800'XL,nws)
    IF nws GT 0 THEN headers(10,ws)=(headers(10,ws) OR 'FFFFF000'XL)+1
    jday=REFORM(headers(0,*))
    hour=REFORM(headers(1,*))
    minute=REFORM(headers(2,*))
    second=REFORM(headers(3,*))
    year=REFORM(headers(4,*))
    glat=REFORM(headers(5,*))/10.
    glon=REFORM(headers(6,*))/10.
    alt=REFORM(headers(7,*))*(6080.*0.3048/1000.)	; naut.miles->km
    clat=REFORM(headers(8,*))/10.
    clon=REFORM(headers(9,*))/10.
    mlat=REFORM(headers(10,*))/10.
    mlon=REFORM(headers(11,*))/10.
    mlt=REFORM(headers(12,*)+(headers(13,*)+headers(14,*)/60.)/60.)

    ; Find the data of interest

    juldate=REPLICATE(julday(1,1,1900+year(0)),nrec)+jday-1
    jd2ymd,juldate,yr,mon,dy
    utsec=(60.*hour+minute)*60.

; These lines applied when stime and etime were in HHMM format.
;   sutsec=(60.*(stime/100)+(stime MOD 100))*60.
;   eutsec=(60.*(etime/100)+(etime MOD 100))*60.

; These lines appy for stime and etime in HHMMSS format.
    stime=ROUND(stime)			; start time, HHMMSS
    shh=stime/10000L
    smm=(stime MOD 10000L)/100L
    sss=stime MOD 100L
    sutsec=(shh*60.+smm)*60.+sss	; start time, sec after 00:00
    fmutsec=60L*FLOOR(sutsec/60)	; even minute before sutsec
    etime=ROUND(etime)			; end time, HHMMSS
    ehh=etime/10000L
    emm=(etime MOD 10000L)/100L
    ess=etime MOD 100L
    eutsec=(ehh*60.+emm)*60.+ess	; end time, sec after 00:00

    nmin=LONG(eutsec)/60L-LONG(sutsec)/60L+1
    IF nmin GT 27 THEN BEGIN
	MESSAGE,STRING(nmin $
		      ,FORMAT='("Selected time interval too long: "' $
			     +',I0," minutes; abort")') $
               ,/INFORMATIONAL
        RETURN
	ENDIF

    wint=WHERE((fmutsec LE utsec) AND (utsec LT eutsec),nwint)

    ; Are there any data of interest?  If so, print them out.

    IF nwint GT 0 THEN BEGIN
	IF KEYWORD_SET(debug) THEN BEGIN
	    fmt=STRING(nwint,FORMAT='("(",I0,"F10.1)")')
	    PRINT,wint,jday(wint),hour(wint),minute(wint),glat(wint) $
		 ,glon(wint),mlat(wint),mlon(wint),mlt(wint) $
		 ,FORMAT=fmt
            ENDIF

	; reopen file and ASSOCiate it to a variable

	OPENR,lun,f(0),/GET_LUN
	file=ASSOC(lun,INTARR(2640))

	; If output desired, open output file

	IF KEYWORD_SET(out) THEN BEGIN
	    oname=name
	    STRPUT,oname,'txt',STRLEN(oname)-3
	    OPENW,olun,oname,/GET_LUN
	    ENDIF

	    ; Open plotting window

	    !ORDER=0
	    gap=8	; pixels between electron and ion panels
	    ; nblocks=9	; old default value

	    ; Load calibration values for this satellite

	    cals=diff(satnum)/1E3	; /keV => /eV
	    esteps=eng(satnum)*1E3	; keV => eV
	    esteps=esteps(0:19)		; same values for e-, ions

	    ; Choose the best UT start and end times, and tick interval

	    getstep,sutsec,eutsec,utstep,psutsec,peutsec ;,/debug
	    sutsec=psutsec
	    eutsec=peutsec
	    nblocks=(eutsec-sutsec)/utstep
	    fact=CEIL((eutsec-sutsec)/600.)

 	    WINDOW,XSIZE=60+600+40+11*!D.X_CH_SIZE $
 		  ,YSIZE=480

 	    dmspct,tcol,bw=bw
 	    whitebot

            datestring=jd2date(juldate(wint(0)),form='d$ n$ Y$')

	    ; Define intermediate arrays to hold plottable data

	    edata=FLTARR(nmin*60,20)
	    idata=edata

	    ; Define and plot colour bar

	    bar=BYTSCL(INDGEN(160),TOP=!D.N_COLORS-1)
	    bar=ROTATE(REBIN(bar,160,20),1)
	    barsize=SIZE(bar)
	    eflux_axis=[10.^(2-(5./14.)),1E7]
            iflux_axis=[10.^(-(5./14.)),1E5]
	    ebarpos=[!D.X_SIZE-barsize(1)-2,!D.Y_SIZE-1*(barsize(2)+gap)-1 $
		    ,!D.X_SIZE-1,!D.Y_SIZE-gap]
	    ibarpos=[!D.X_SIZE-barsize(1)-2,!D.Y_SIZE-2*(barsize(2)+gap)-1 $
		    ,!D.X_SIZE-1,!D.Y_SIZE-1*(barsize(2)+gap)-gap+1]
 	    PLOT_IO,[0,1],eflux_axis,/DEVICE,/NODATA,POSITION=ebarpos $
 		   ,XRANGE=[0,1],XSTYLE=1,XTICKLEN=0.0,XTICKS=1 $
 		   ,XMINOR=0,XTICKV=INDGEN(2),XTICKNAME=REPLICATE(" ",2) $
 		   ,YRANGE=eflux_axis,YSTYLE=1,YTICKLEN=0.0,YTICKS=2 $
 		   ,YMINOR=1,YTICKV=[1E2,1E4,1E6],YTICKFORMAT='power10' $
 		   ,YTITLE='!17!CELECTRON FLUX',YCHARSIZE=1.2 $
 		   ,COLOR=tcol
 	    PLOT_IO,[0,1],iflux_axis,/DEVICE,/NODATA,POSITION=ibarpos $
 		   ,XRANGE=[0,1],XSTYLE=1,XTICKLEN=0.0,XTICKS=1 $
 		   ,XMINOR=0,XTICKV=INDGEN(2),XTICKNAME=REPLICATE(" ",2) $
 		   ,YRANGE=iflux_axis,YSTYLE=1,YTICKLEN=0.0,YTICKS=2 $
 		   ,YMINOR=1,YTICKV=[1E0,1E2,1E4],YTICKFORMAT='power10' $
 		   ,YTITLE='!CION FLUX',YCHARSIZE=1.2 $
 		   ,COLOR=tcol,/NOERASE
 	    TV,bar,!D.X_SIZE-barsize(1)-1,!D.Y_SIZE-1*(barsize(2)+gap),/DEVICE
 	    TV,bar,!D.X_SIZE-barsize(1)-1,!D.Y_SIZE-2*(barsize(2)+gap),/DEVICE

	    ; Set up plot frames

	    taxis=sutsec+FINDGEN(eutsec-sutsec+1)
	    eaxis=[(esteps(19)^1.5/esteps(18)^0.5) $
		  ,(esteps(0)^1.5/esteps(1)^0.5)]

	    ; Note that the time axis length is fixed at 600 pixels

	    epos=[60-1,!D.Y_SIZE-1*(barsize(2)+gap)-1 $
                 ,60+600+1,!D.Y_SIZE-gap]
 	    PLOT_IO,taxis,eaxis,/DEVICE,/NODATA,POSITION=epos $
 		   ,XRANGE=[sutsec,eutsec],XSTYLE=1,XTICKLEN=-0.02 $
 		   ,XTICKS=nblocks,XMINOR=0 $
 		   ,XTICKV=FINDGEN(nblocks+1)*utstep+sutsec $
 		   ,XTICKNAME=REPLICATE(" ",nblocks+1) $
 		   ,YRANGE=eaxis,YSTYLE=1,YTICKLEN=0.0 $
 		   ,YTICKS=2,YMINOR=0,YTICKV=[1E2,1E3,1E4] $
 		   ,YTICKFORMAT='power10',YTITLE='!CENERGY (eV)' $
 		   ,YCHARSIZE=1.2,COLOR=tcol,/NOERASE

	    ; Note again that the time axis length is fixed at 600 pixels

	    ipos=[60-1,!D.Y_SIZE-2*(barsize(2)+gap)-1 $
		 ,60+600+1,!D.Y_SIZE-1*(barsize(2)+gap)-gap+1]
 	    PLOT_IO,taxis,eaxis,/DEVICE,/NODATA,POSITION=ipos $
 		   ,XRANGE=[sutsec,eutsec],XSTYLE=1,XTICKLEN=-0.02 $
 		   ,XTICKS=nblocks,XMINOR=0 $
 		   ,XTICKV=FINDGEN(nblocks+1)*utstep+sutsec $
 		   ,XTICKNAME=REPLICATE(" ",nblocks+1) $
 		   ,YRANGE=eaxis,YSTYLE=1,YTICKLEN=0.0 $
 		   ,YTICKS=2,YMINOR=0,YTICKV=[1E2,1E3,1E4] $
 		   ,YTICKFORMAT='power10',YTITLE='!C'+datestring $
 		   ,YCHARSIZE=1.2,COLOR=tcol,/NOERASE

            satlabel='F'+STRTRIM(satnum,2)
 	    XYOUTS,1,!D.Y_SIZE-1-!D.Y_CH_SIZE,satlabel $
 		  ,ALIGNMENT=0,CHARSIZE=1.2,/DEVICE,COLOR=tcol
		  
	; Set up the unscrambling array for the count rate codes
	resort=ROTATE(REFORM(INDGEN(20),4,5),5)

	IF NOT KEYWORD_SET(out) THEN BEGIN
            resp=''
            READ,'Do you want to print to the screen? (Y/[N]) ' $
		,resp
            IF STRUPCASE(resp) EQ 'Y' THEN iocmd='PRINT'
            ENDIF ELSE iocmd='PRINTF,olun'

	FOR i=0,nwint-1 DO BEGIN
	    IF KEYWORD_SET(debug) THEN PRINT,wint(i)
	    IF KEYWORD_SET(out) THEN PRINTF,olun,wint(i)
	    nums=file(wint(i))
	    IF !VERSION.OS EQ 'RISC/os' THEN BYTEORDER,nums

	    data=REFORM(nums(15:2594),43,60)
	    times=data(0:2,*)
	    codes=data(3:42,*)
	    cts=2^(codes/32)*((codes MOD 32)+32)-33

	    ; Flag negative data

	    IF MIN(cts) LT 0 THEN BEGIN
		wneg=WHERE(cts LT 0,nwneg)
		one2two,wneg,cts,xe,yt
		cts(wneg)=0
		ENDIF

	    ; Separate "electron" and "ion" data and unscramble them

	    e_cts=cts(resort,*)
	    i_cts=cts(20+resort,*)

	    ; Convert counts to fluxes

	    e_flux=e_cts*REBIN(cals(0:19),20,60)
	    i_flux=i_cts*REBIN(cals(20:39),20,60)

	    ; Save 1 minute of fluxes

	    edata(utsec(wint(i))-fmutsec,0)=TRANSPOSE(e_flux)
	    idata(utsec(wint(i))-fmutsec,0)=TRANSPOSE(i_flux)

	    IF N_ELEMENTS(iocmd) EQ 1 THEN BEGIN
 	        result=EXECUTE(iocmd+",nums(0:14),FORMAT='(15I6)'")
 	        result=EXECUTE(iocmd+",e_cts,FORMAT='(20I6)'")
 	        result=EXECUTE(iocmd+",FORMAT='()'")
 	        result=EXECUTE(iocmd+",i_cts,FORMAT='(20I6)'")
 	        result=EXECUTE(iocmd+",FORMAT='()'")
		ENDIF

	    ENDFOR

	    FREE_LUN,lun

	; Truncate edata and idata arrays to contain only data
	; between sutsec and eutsec

	allut=makex(fmutsec,60L*CEIL(eutsec/60)-1,1)
	toplot=WHERE((sutsec LE allut) AND (allut LT eutsec),ntoplot)
	allut=allut(toplot)
	edata=edata(toplot,*)
	idata=idata(toplot,*)

	; Average data if more than 9 minutes are to be presented.

            IF fact GT 1 THEN BEGIN
	        pedata=REBIN(edata,ntoplot/fact,20)
	        pidata=REBIN(idata,ntoplot/fact,20)
	        ENDIF ELSE BEGIN
		pedata=edata
		pidata=idata
		ENDELSE

	    ; Convert from fluxes to bytes for plotting

	    ebytes=BYTSCL(ALOG10(pedata>10) $
			 ,MIN=1.6429,MAX=7.,TOP=!D.N_COLORS-1)
	    ibytes=BYTSCL(ALOG10(pidata>0.1) $
			 ,MIN=-0.3571,MAX=5.,TOP=!D.N_COLORS-1)

	    ; Sample array to fill plot space

	    tsamp=FINDGEN(600)*((ntoplot/fact)/599.)#REPLICATE(1,1,160)

	    ; This block performs correct sampling over the plotted
	    ; energy range

	    emid=[eaxis(1),SQRT(esteps(0:18)*esteps(1:19)),eaxis(0)]
	    es=ALOG10(emid)		; energy step
	    cn=FINDGEN(21)-0.5		; channel number
	    esi=(FINDGEN(160)/159.)*ALOG10(eaxis(1)/eaxis(0))+ALOG10(eaxis(0))
	    cni=(ROUND(interpol(cn,es,esi))>0)<19
	    esamp=REPLICATE(1,600)#cni
	    eout=ebytes(tsamp,esamp)
	    iout=ibytes(tsamp,esamp)

	    ; Plot the data

 	    TV,eout,60,!D.Y_SIZE-(barsize(2)+gap),/DEVICE
 	    TV,iout,60,!D.Y_SIZE-2*(barsize(2)+gap),/DEVICE

	    ; Get extra (and more accurate) position data using an
	    ; external orbital calculating program and accurate
	    ; orbital elements.

	    ; First define the inputs to the PLOTORB routine

	    satname='DMSP B5D2-'+STRING(satnum-5,FORMAT='(I1)')
	    podate=(((yr MOD 100)*100L+mon)*100L+dy)(wint(0))
	    sechms,sutsec,shh,smm,sss
	    sechms,eutsec,ehh,emm,ess
	    IF ess EQ 0 THEN ess=1	; kludge to make 'orbit' work
	    postart=(shh*100L+smm)*100L+sss
	    poendt=(ehh*100L+emm)*100L+ess
	    point=utstep/60.
	    potype=2

	    ; Call PLOTORB to generate footprint positions at the 
	    ; default altitude of 250 km

            potype=4
	    IF N_ELEMENTS(fpalt) EQ 0 THEN fpalt=250.
	    plotorb,satname,podate,postart,poendt,point $
		   ,fptime,fplat,fplong,fpht $
		   ,type=potype,/noplot,fpalt=fpalt
	    fptime=ROUND(fptime)
            fplong=(fplong+360.0) MOD 360

	    ; Convert these positions to PACE geomagnetic
	    ; coordinates with the MCONVERT routine; first
	    ; set up the inputs

	    nfp=N_ELEMENTS(fplat)
	    myr=podate/10000L
	    mmo=(podate MOD 10000L)/100
	    mdy=(podate MOD 100)
	    mhr=fptime/10000L
	    mmi=(fptime MOD 10000L)/100
	    msc=(fptime MOD 100)
	    pgmlat=FLTARR(nfp)
	    pgmlon=pgmlat
	    pgmlt=pgmlat

	    ; Then call MCONVERT for each position (it's in
	    ; FORTRAN and can't handle vector data!)

	    IF nfp NE nblocks+1 THEN BEGIN
		MESSAGE,'Wrong number of PACE positions calculated!' $
		       ,/INFORMATIONAL
		PRINT,nblocks,nfp
		STOP
		ENDIF

            FOR ipt=0,nfp-1 DO BEGIN
		mconvert,myr,mmo,mdy,mhr(ipt),mmi(ipt),msc(ipt) $
			,fplat(ipt),fplong(ipt),fpht(ipt),1 $
			,tmlat,tmlon,tmlt
                pgmlat(ipt)=tmlat
		pgmlon(ipt)=(tmlon+360.) MOD 360
		pgmlt(ipt)=tmlt
		ENDFOR

	    ; Decorate with times and places

	    pout=sutsec+FINDGEN(nfp)*60*point
	    ydev=!D.Y_SIZE-2*(barsize(2)+gap)-!D.Y_CH_SIZE
	    IF (sss EQ 0) AND (utstep GE 60) $	; fix no. of HHMMSS
		THEN hmschar=5 $		; characters displayed:
		ELSE hmschar=8			; HH:MM or HH:MM:SS

	    FOR i=0,nfp-1 DO BEGIN
 	        labelpos=CONVERT_COORD(pout(i),eaxis(0) $
 				      ,/DATA,/TO_DEVICE)
                xdev=labelpos(0)

;		IF i EQ 0 THEN PRINT,pout,pgmlat,pgmlt,fplat,fplong

 	        XYOUTS,xdev,ydev-1.5*!D.Y_CH_SIZE $
 		      ,STRING(LONG(pout(i)),FORMAT='(I5)') $
 		      ,/DEVICE,ALIGNMENT=0.5,COLOR=tcol
 	        XYOUTS,xdev,ydev-3.0*!D.Y_CH_SIZE $
 		      ,STRMID(strsec(pout(i),/hours),0,hmschar) $
 		      ,/DEVICE,ALIGNMENT=0.5,COLOR=tcol
 	        XYOUTS,xdev,ydev-4.5*!D.Y_CH_SIZE $
 		      ,STRING(pgmlat(i),FORMAT='(F6.2)') $
 		      ,/DEVICE,ALIGNMENT=0.5,COLOR=tcol
 	        XYOUTS,xdev,ydev-6.0*!D.Y_CH_SIZE $
 		      ,STRING(pgmlt(i),FORMAT='(F6.2)') $
 		      ,/DEVICE,ALIGNMENT=0.5,COLOR=tcol
 	        XYOUTS,xdev,ydev-7.5*!D.Y_CH_SIZE $
 		      ,STRING(fplat(i),FORMAT='(F6.2)') $
 		      ,/DEVICE,ALIGNMENT=0.5,COLOR=tcol
 	        XYOUTS,xdev,ydev-9.0*!D.Y_CH_SIZE $
 		      ,STRING(fplong(i),FORMAT='(F6.2)') $
 		      ,/DEVICE,ALIGNMENT=0.5,COLOR=tcol
                ENDFOR

 	    keypos=CONVERT_COORD(eutsec,eaxis(0),/DATA,/TO_DEVICE)
 	    xdev=keypos(0)+6.5*!D.X_CH_SIZE

 	    XYOUTS,xdev,ydev-1.5*!D.Y_CH_SIZE,'UT(SEC)' $
 	          ,/DEVICE,ALIGNMENT=0,COLOR=tcol
	    IF (sss EQ 0) AND (utstep GE 60) $
		THEN hhstr='HH:MM' $
		ELSE hhstr='HH:MM:SS'
 	    XYOUTS,xdev,ydev-3.0*!D.Y_CH_SIZE,hhstr $
 	          ,/DEVICE,ALIGNMENT=0,COLOR=tcol
 	    XYOUTS,xdev,ydev-4.5*!D.Y_CH_SIZE,'PGMLAT' $
 	          ,/DEVICE,ALIGNMENT=0,COLOR=tcol
 	    XYOUTS,xdev,ydev-6.0*!D.Y_CH_SIZE,'PGMLT' $
 	          ,/DEVICE,ALIGNMENT=0,COLOR=tcol
 	    XYOUTS,xdev,ydev-7.5*!D.Y_CH_SIZE,'FPLAT' $
 	          ,/DEVICE,ALIGNMENT=0,COLOR=tcol
 	    XYOUTS,xdev,ydev-9.0*!D.Y_CH_SIZE,'FPLON' $
 	          ,/DEVICE,ALIGNMENT=0,COLOR=tcol

	; Do we want to plot spectra?

	IF KEYWORD_SET(spectra) THEN BEGIN

	    psout=''
	    READ,'Plot spectra to a PostScript file? (Y/[N]) ',psout
	    IF STRUPCASE(psout) EQ 'Q' THEN RETURN
	    ps=(STRUPCASE(psout) EQ 'Y')

	    PRINT,'Click on the species and time for which spectra' $
		 +' are desired:'
            CURSOR,xtime,yspecies,/DEVICE,/UP

 	    REPEAT BEGIN
 	        okay= ((epos(0) LT xtime) AND (xtime LT epos(2))) AND $
 		     (((ipos(1) LT yspecies) AND (yspecies LT ipos(3))) OR $
 		      ((epos(1) LT yspecies) AND (yspecies LT epos(3))))
 		IF NOT okay THEN BEGIN
 		    MESSAGE,'Cursor out of range, try again!' $
 			   ,/INFORMATIONAL
                    CURSOR,xtime,yspecies,/DEVICE,/UP
 		    ENDIF
                ENDREP UNTIL okay

            electron=(epos(1) LT yspecies) AND (yspecies LT epos(3))

;	    READ,'Enter 1 for e- or 0 for ions: ',electron

	    IF electron THEN BEGIN
		spdata=edata 
		spcals=cals(0:19)
		species='!17Electron '
		MESSAGE,'Plotting electron spectra',/INFORMATIONAL
		ENDIF ELSE BEGIN
		spdata=idata
		spcals=cals(20:39)
		species='!17Ion '
		MESSAGE,'Plotting ion spectra',/INFORMATIONAL
		ENDELSE

 	    sput=ROUND((FLOAT(xtime-epos(0))/(epos(2)-epos(0))) $
 		      *(eutsec-sutsec)+sutsec)

;	    PRINT,sutsec,eutsec $
;		 ,FORMAT='("Data span ",F6.0," - ",F6.0," s;")'

 	    fmt='("You have selected ",F6.0," s (",A,") UT.")'
 	    PRINT,STRING(sput,strsec(sput),FORMAT=fmt)
	    READ,'Enter the exact center time and bin width (UT s): ' $
		,sput,binwid
	    nplots=9
	    tol=(nplots/2.)*binwid
	    spnum=WHERE(ABS(allut-sput) LT tol,nspnum)

	    IF ps THEN BEGIN
		spstart=allut(spnum(0))
		psspname=STRMID(f(0),0,STRLEN(f(0))-4)+'_'
		IF electron THEN spchar='e' ELSE spchar='i'
		psspname=psspname+spchar+'_'
		sechms,spstart,h,m,s
		sphms=(h*100L+m)*100L+s
		psspname=psspname+STRING(sphms,binwid $
					,FORMAT='(I6.6,"_",I0,".ps")')
                PRINT,psspname
		psopen,psspname,/long
		ENDIF

	    IF nspnum EQ 0 THEN BEGIN
		MESSAGE,'No spectra found; abort',/INFORMATIONAL
		RETURN
		ENDIF ELSE BEGIN

		spect=FLTARR(9,20)
		spj=FLTARR(9)
		spje=spj
		spavge=spj
		sptime=STRARR(9)
		bumspect=BYTARR(9)
		FOR i=0,(nspnum-1)/binwid DO BEGIN
		    xp=INDGEN(binwid)+i*binwid
		    inspnum=WHERE((0 LE xp) AND (xp LT nspnum),ninspnum)
		    IF ninspnum GT 0 THEN BEGIN
			xp=xp(inspnum)
		        gxp=WHERE(TOTAL(spdata(spnum(xp),*),2) GT 0,ngxp)
		        IF ngxp GT 0 THEN BEGIN
			    ; gxp=gxp+i*binwid
			    xp=xp(gxp)

;			    k=GET_KBRD(1)

		            spect(i,*)=TOTAL(spdata(spnum(xp),*),1)/ngxp
			    spint,esteps,REFORM(spect(i,*)),je,avge,j
			    spj(i)=j
			    spje(i)=je
			    spavge(i)=avge/1E3		; eV => keV

;			    k=GET_KBRD(1)

			    ts=STRMID(strsec(allut(spnum(xp(0))),/hours),0,8)
			    ts=ts+' - '
			    ts=ts+STRMID(strsec(allut(spnum(xp(ngxp-1))) $
					       ,/hours) $
					,0,8)
			    ts=ts+' UT'
			    sptime(i)=ts
;			    IF ngxp LT binwid THEN STOP
			    ENDIF			; (ngxp GT 0)
			ENDIF ELSE bumspect(i)=1	; (ninspnum GT 0)
		    ENDFOR				; (i=...)

		gspect=WHERE(TOTAL(spect,2) GT 0,nspect)
		maxj=MAX(spect(WHERE(spect GT 0)),MIN=minj)
		wbum=WHERE(bumspect,nwbum)
                IF nwbum GT 0 THEN $
		    FOR i=0,nwbum-1 DO spect(wbum(i),*)=REPLICATE(minj,20)

                IF !D.NAME EQ 'X' THEN WINDOW,2,XSIZE=800,YSIZE=1000
		multiplot,[3,3],/init
		!X.RANGE=[10.,9E4]
		!X.STYLE=1
		loj=10.^FLOOR(ALOG10(minj))
		hij=10.^CEIL(ALOG10(maxj))
		!Y.RANGE=[loj,hij]
		!Y.STYLE=1
		ptit=[' ' $
		     ,'!17DMSP '+satlabel+' Precipitating '+species $
		     +'Spectra, '+datestring]
		xtit=  [' ',species+'Energy (eV)']
		xtickf=['nullst','power10']
		ytit=  [' ','!17dJ/dE (cm!e2!n s sr eV)!e-1!n']
		ytickf=['nullst','power10']

		hiE=INDGEN(10)		; plot two halves of spectra
		loE=hiE+10		; separately 

		FOR i=0,nspect-1 DO BEGIN
		    multiplot
		    ppt=(i EQ 1)
		    pxt=(i/3 EQ 2)
		    pyt=(i MOD 3) EQ 0
		    PLOT_OO,esteps(hiE),spect(gspect(i),hiE) $
			   ,PSYM=-4,COLOR=2,TITLE=ptit(ppt) $
			   ,XTITLE=xtit(pxt) $
			   ,XTICKFORMAT=xtickf(pxt) $
			   ,YTITLE=ytit(pyt) $
			   ,YTICKFORMAT=ytickf(pyt)
		    OPLOT,esteps(loE),spect(gspect(i),loE) $
			 ,PSYM=-4,COLOR=2
                    OPLOT,esteps,spcals,COLOR=2
		    datacorn=[TRANSPOSE(10.^!X.CRANGE) $
			     ,TRANSPOSE(10.^!Y.CRANGE) $
			     ,TRANSPOSE([0.0,0.0])]
		    devcorn=FLOAT(CONVERT_COORD(datacorn $
					       ,/DATA,/TO_DEVICE))
                    nchars=(devcorn(0,1)-devcorn(0,0))/!D.X_CH_SIZE
		    !P.CHARSIZE=(nchars/22)<1
		    XYOUTS,devcorn(0,1)-!D.X_CH_SIZE $
			  ,devcorn(1,1)-1.5*!D.Y_CH_SIZE $
			  ,sptime(gspect(i)) $
			  ,/DEVICE,ALIGNMENT=1,COLOR=2
;		    IF spavge(i) LT 9.99 $
;			THEN aefmt='("<E>=",F4.2," keV")' $
;			ELSE aefmt='("<E>=",F4.1," keV")'
;		    XYOUTS,devcorn(0,0)+!D.X_CH_SIZE $
;			  ,devcorn(1,0)+1.75*!D.Y_CH_SIZE $
;			  ,STRING(spavge(gspect(i)),FORMAT=aefmt) $
;			  ,/DEVICE,ALIGNMENT=0,COLOR=2
		    CASE 1 OF
                       (spje(i) LT 0.001): jefmt='(F6.4'
(0.001 LE spje(i)) AND (spje(i) LE 9.999): jefmt='(F5.3'
                                     ELSE: jefmt='(F6.3'
			ENDCASE
		    jefmt=jefmt+'," erg cm!e-2!n s!e-1!n sr!e-1!n")'
		    XYOUTS,devcorn(0,0)+!D.X_CH_SIZE $
			  ,devcorn(1,0)+0.5*!D.Y_CH_SIZE $
			  ,STRING(spje(gspect(i)),FORMAT=jefmt) $
			  ,/DEVICE,ALIGNMENT=0,COLOR=2
                    ENDFOR

		IF ps THEN psclose

                ENDELSE

	    ENDIF

        ENDIF

    ; That's it for now

    ENDELSE

    IF N_ELEMENTS(iocmd) EQ 1 THEN BEGIN
        result=EXECUTE(iocmd+",ALOG10(TRANSPOSE(edata>10)),FORMAT='(20F6.2)'")
        result=EXECUTE(iocmd+",FORMAT='()'")
        result=EXECUTE(iocmd+",ALOG10(TRANSPOSE(idata>0.1)),FORMAT='(20F6.2)'")
        result=EXECUTE(iocmd+",FORMAT='()'")
	ENDIF
    IF KEYWORD_SET(out) THEN FREE_LUN,olun

RETURN
END