pro DISTRIB, vplt = IVplt, tplt = ITYplt

@compile_opt.pro			; On error, return to caller
common	DISARR	, ISC, DVEL, DQ, Points, IVsave, ITYsave
common  MAPDIR	, direc

if !d.NAME eq !ThePrinter then goto, PLOT_HISTOGRAMS

if IsType(IVplt,/def)*IsType(ITYplt,/def) then goto, PLOT_HISTOGRAMS

InitVar, direc, 'ud1:[plh.data.hb.square]'
echo, 'Directory ?', direc
cd, direc  &  dir

; Interactive prompts
; The data file names are of type A0MAP003_5.0300
;	A	= spacecraft
;	0	= data type nV^0
;	003	= start Carrington rotation
;	_5	= # rotations averaged
;	.0300	= reference distance 0.3 AU

ISC = 1  &  echo, 'Spacecraft ( 1=A , 2=B )', ISC, 1,2
MAP_MASK = 'A0MAP003_5.0300'		; Sample Helios map file
strput, MAP_MASK, strmid('0AB',ISC,1), 0

; Set bin width for solar wind velocity (defines # velocity bins)

DVEL   = 1000					; velocity bin width (km/s)
print  &  echo, 'Velocity bin width (km/s)', DVEL, 10., 1000.

NVEL   = 1+round(1000./DVEL)			; # velocity classes
echo, '# velocity bins : '+strcompress(NVEL,/remove_all)

NTYP   = 4					; # data types
Qmax   = [70., 240., 120., 300.]		; Max. flux values
NQ     = 20					; # flux bins
DQ     = Qmax/NQ				; Flux bin width

; Each data file is identified by a start Carrington rotation and the 
; # rotations summed together.

echo, 'Specify maps to be summed'
ICA = 3  &  ISU = 5  &  NCA = -1  &  NSU = -1  &  MAPS = '0'

echo, 'Carr. rotation (0=End selection)', ICA, 0, 999
while ICA ne 0 do begin
    echo, 'Total # Carrington rotations', ISU, 1,9
    strput, MAP_MASK, string(format='(I3.3)',ICA)  , 5
    strput, MAP_MASK, strcompress(ISU,/remove_all) , 9
    MAP_ = file_search( MAP_MASK )
    if MAP_[0] ne '' then begin
	NCA = [NCA,ICA]  &  NSU = [NSU,ISU]  &  MAPS = [MAPS,MAP_[0]]
	message, /info, 'File '+MAP_[0]
    endif else message, /info, 'File '+MAP_MASK+' not found'
    echo, 'Carr. rotation (0=End selection)', ICA, 0, 999
endwhile

if n_elements(NCA) eq 1 then message,'No maps specified'
NCA = NCA[1:*]  &  NSU = NSU[1:*]  &  MAPS = MAPS[1:*]
CTY = strpos ( MAPS , ']' ) + 2

; Read the start times for all HELIOS Carrington rotations into array JDCAR

NCAR = 151  &  JDCAR = dblarr(NCAR)  &  I = 0  &  JD = 0.d

openr, /get_lun, IU, filepath(root='$AUX','HELIOS'+strcompress(ISC,/remove_all)+'.CAR')
while not eof(IU) do begin
    readf, IU, format='(I7,F20.9)', J, JD  &  I = I+1  &  JDCAR[I] = JD
endwhile
free_lun, IU

; Create arrays to contain binned data

Points = intarr(NVEL,NQ,NTYP)	; # points in each velocity class
help, points

NMAP = n_elements(NCA)		; # Helios maps
VEL_FILE = 'V1612_8.030'	; Sample velocity file name
NX = 37  &  NY = 19  &  NYH = 10  &  JSC = (ISC-1)*(NYH-1)

; Weight factor is cos (latitude)

DY = (!pi/2)/(NYH-1)

for IMA=0,NMAP-1 do begin		; Loop over all selected maps

; Locate the appropriate SANIPS velocity file

    JD_EAR = 0.5*( JDCAR[NCA[IMA]+NSU[IMA]-1] + JDCAR[NCA[IMA]] )
    JULIAN, 1, IYR_EAR, DOY_EAR, JD_EAR
    N_EAR = N_CARRINGTON(IYR_EAR,DOY_EAR)
print, iyr_ear, doy_ear, n_ear
    strput, VEL_FILE, string ( format='(I4.4)', N_EAR ) , 1
    CFOUND = (file_search( filepath(root='$IPS',strmid(VEL_FILE,0,5)+'*.030' )))[0]

; If the first guess was wrong search for adjacent Carrington rotations

    I = 0
    while CFOUND eq '' and I lt 7 do begin
   	I = I+1
   	strput, VEL_FILE, string ( format='(I4.4)', N_EAR-I ) , 1
   	CFOUND = (file_search( filepath(root='$IPS',strmid(VEL_FILE,0,5)+'*.030' )))[0]
   	if CFOUND eq '' then begin
	    strput, VEL_FILE, string ( format='(I4.4)', N_EAR+I ) , 1
	    CFOUND = (file_search( filepath(root='$IPS',strmid(VEL_FILE,0,5)+'*.030' )))[0]
	endif
    endwhile
    if CFOUND eq '' then message, 'no matching IPS velocity file >> '+VEL_FILE
    message, /info, 'IPS velocity file >> '+CFOUND

; Read velocity file. Extract hemisphere corresponding to ISC.
; For Helios 1 (ISC=1) flip array upside down (so that IY=0 becomes equator)

    stat = flt_read ( CFOUND, VMAP , nx=NY , ny=NX, /stop )
    VMAP = transpose ( VMAP )
    VMAP = VMAP[*,JSC:JSC+NYH-1]
    if ISC eq 1 then VMAP = rotate ( VMAP , 7 )		; Y --> -Y

; Velocity range in class IV is (IV-1)*DVEL to IV*DVEL, (IV=1,NVEL)

    W = where ( VMAP gt 0 , Cnt )
    if Cnt gt 0 then VMAP[W] = 1+fix ( VMAP[W]/DVEL )

    for ITY=0,NTYP-1 do begin				; Loop over data types

; Insert data type into file name, open file

	MAP_ = MAPS[IMA]
	strput, MAP_ , string ( format='(I1.1)' , ITY ), CTY[IMA]
	MAP_ = file_search ( MAP_ )  &  MAP_ = MAP_[0]
	if MAP_ eq '' then message, 'Helios map file not found'

; Read Helios map file
; For Helios 1 (ISC=1) flip array upside down (so that IY=0 becomes equator)

	stat = flt_read ( MAP_ , QMAP, nx=NYH , ny=2*NX , /stop )
	QMAP = transpose ( QMAP )
	QMAP = QMAP[0:NX-1,*]
	if ISC eq 1 then QMAP = rotate ( QMAP , 7 )

; Assign flux bin values
	W = where ( QMAP gt -1000. , Cnt )
	if Cnt gt 0 then QMAP[W] = fix ( QMAP[W]/DQ[ITY] )

; Set invalid flux bins ( value -1000. to -1)

	QMAP = QMAP > (-1)

; Register invalid bins (no velocity or no flux) in velocity class 0

	Invalid = where ( VMAP eq 0 or QMAP lt 0 , Cnt )
	if Cnt gt 0 then Points [0,0,ITY] = Points[0,0,ITY] + Cnt

; Calculate # points and sum of weights in each velocity class
; Velocity range in class IV is (IV-1)*DVEL to IV*DVEL

	for IV=1,NVEL-1 do begin		; Loop over velocity classes
	    for IQ = 0,NQ-1 do begin		; Loop over flux classes
		Valid = where ( VMAP eq IV and QMAP eq IQ , Cnt )
		if Cnt gt 0 then Points [IV,IQ,ITY] = Points [IV,IQ,ITY] + Cnt
	    endfor
	endfor

    endfor

endfor

IVplt = 1  &  ITYplt = 0

PLOT_HISTOGRAMS:

if !d.NAME eq !ThePrinter then begin IVplt = IVsave  &  ITYplt = ITYsave  &  endif
if IVplt eq 0 then ITYplt = 0

if IVplt ne 0 then begin
    MTITLE	=	strcompress( (IVplt-1)*round(DVEL), /remove_all ) + $
			' - ' + strcompress( IVplt*round(DVEL), /remove_all ) + $
			' km/s'
    YTITLE	=	'TYPE '+strcompress(ITYplt)
endif else begin
    MTITLE	=	'invalid bins'
    YTITLE	=	'invalid bins'
endelse

if max( Points[IVplt,*,ITYplt] ) eq 0 then		$
    message, /info, 'No bins available'+		$
	', IVplt  = '+strcompress(IVplt , /remove_all) +$
	', ITYplt = '+strcompress(ITYplt, /remove_all)	$
else begin
    bargraph, Points[IVplt,*,ITYplt] , 0., DQ[ITYplt],	$
	mtitle=MTITLE, ytitle=YTITLE
    xyouts, /normal, .4,0., systime(), charsize=.75	; Current date/time
endelse

setup_ooze

;davis, 1  &  twin, 2

IVsave = IVplt  &  ITYsave = ITYplt

return  &  end