;
pro rdprotons,ps=ps
;
; Read ascii data file with Jackman proton data from 5 events:
;
;rdfile = 'jackman_5spe.dat'
rdfile = 'jackman.jan05.dat'
openr,lurd,rdfile,/get_lun
header = ''
;
; Determine number of events:
;
nevents = 0
while not eof(lurd) do begin
  readf,lurd,header
  if strmid(header,0,1) eq '>' then begin
    nevents = nevents+1
  endif
endwhile
print,'nevents=',nevents
ntimes = intarr(nevents) ; number of times in each event
point_lun,lurd,0
;
; Find number of times in each event:
;
ievent = -1
while not eof(lurd) do begin
  readf,lurd,header
  if strmid(header,0,1) eq '>' then begin ; new event
    for i=1,8 do begin
      readf,lurd,header
;     print,header
    endfor
    ievent = ievent+1
    ntimes[ievent] = 0
  endif else begin
    ntimes[ievent] = ntimes[ievent]+1
  endelse
endwhile
;
; Allocate data:
;
print,'ntimes = ',ntimes
mxtimes = max(ntimes)
fe = fltarr(mxtimes,nevents)
eo = fltarr(mxtimes,nevents)
iday = intarr(mxtimes,nevents)
hour = fltarr(mxtimes,nevents)
;
; Read data:
;
hdrdates = strarr(nevents)
ievent = -1
point_lun,lurd,0
fe_rd=0. & eo_rd=0. & iday_rd=0 & hour_rd=0.
while not eof(lurd) do begin
  readf,lurd,header
  if strmid(header,0,1) eq '>' then begin ; new event
    for i=1,8 do begin
      readf,lurd,header
;     print,header
      if i eq 1 then hdrdates[ievent+1] = header
    endfor
    ievent = ievent+1
    itime = -1
    ntime = ntimes[ievent]
    for i=0,ntime-1 do begin
      readf,lurd,fe_rd,eo_rd,iday_rd,hour_rd
      itime = itime+1
      fe[itime,ievent] = fe_rd
      eo[itime,ievent] = eo_rd
      iday[itime,ievent] = iday_rd
      hour[itime,ievent] = hour_rd
    endfor
    print,hdrdates[ievent]
;   print,'Fe = ' & print,fe[0:ntime-1,ievent]
;   print,'eo = ' & print,eo[0:ntime-1,ievent]
;   print,'iday = ' & print,iday[0:ntime-1,ievent]
;   print,'hour = ' & print,hour[0:ntime-1,ievent]
  endif ; new event
endwhile
free_lun,lurd
;
; Plots:
;
if keyword_set(ps) then set_plot,'ps'
!p.multi = [0,2,3,1,0]
for i=0,nevents-1 do begin
  ntime = ntimes[i]
  title = strmid(hdrdates[i],27)
  dday = iday[0:ntime-1,i]+hour[0:ntime-1,i]/60.
  plot,dday,fe[0:ntime-1,i],title=title,$
    ytitle='Fe energy flux (MeV/cm2/s)',$
    xtitle='Day',xrange=[dday[i],dday[ntime-1]]
endfor
!p.multi = [0,2,3,1,0]
for i=0,nevents-1 do begin
  ntime = ntimes[i]
  title = strmid(hdrdates[i],27)
  dday = iday[0:ntime-1,i]+hour[0:ntime-1,i]/60.
  plot,dday,eo[0:ntime-1,0],title=title,$
    ytitle='Eo characteristic energy (MeV)',$
    xtitle='Day',xrange=[dday[0],dday[ntime-1]]
endfor
if keyword_set(ps) then device,/close
;
; Write f90 source file with array constructors
; (this can be refined by hand later for import to the model).
;
filew = 'protons.F'
openw,luw,filew,/get_lun
printf,luw,'!'
printf,luw,'      subroutine protons'
printf,luw,'      implicit none'
printf,luw,format="('      integer,parameter,dimension(',i1,') :: ntimes = (/',5(i3,','),'/)')",nevents,ntimes
printf,luw,format="('      real,dimension(',i3,',',i1,') :: fe,eo,hour')",$
  mxtimes,nevents
printf,luw,format="('      integer,dimension(',i3,',',i1,') :: iday')",$
  mxtimes,nevents
;
for i=0,nevents-1 do begin
  ntime = ntimes[i]
  printf,luw,'!'
  printf,luw,format="('! ',a)",hdrdates[i] 
;
; Array constructor for fe:
  printf,luw,format="('      fe(1:',i3,',',i1,')=(/')",ntimes[i],i+1
  printf,luw,format="(('     |',6(e10.3,',')))",fe[0:ntime-1,i]
  printf,luw,'     | /)'
;
; Array constructor for eo:
  printf,luw,format="('      eo(1:',i3,',',i1,')=(/')",ntimes[i],i+1
  printf,luw,format="(('     |',6(e10.3,',')))",eo[0:ntime-1,i]
  printf,luw,'     | /)'
;
; Array constructor for iday:
  printf,luw,format="('      iday(1:',i3,',',i1,')=(/')",ntimes[i],i+1
  printf,luw,format="(('     | ',12(i4,',')))",iday[0:ntime-1,i]
  printf,luw,'     | /)'
;
; Array constructor for hour:
  printf,luw,format="('      hour(1:',i3,',',i1,')=(/')",ntimes[i],i+1
  printf,luw,format="(('     | ',12(f4.1,',')))",hour[0:ntime-1,i]
  printf,luw,'     | /)'
endfor
printf,luw,'      end subroutine protons'
close,luw
end ; rdprotons