pro read_mrr_raw_txt


U=''
Frequency=replicate(!values.f_nan,[64,32])
dummy=replicate(!values.F_NAN,[64,32])
first=-1

;  Get the data from the txt file
file='1211.raw.txt'
openr,lun,file,/get_lun

;while not eof(lun) do begin
for ii=0,(4*67) do begin
  line=''
  readf,lun,line
  ;print,line
  indicator=''
  reads,line,indicator,format='(A1)'
  ;*** M *** TIME
  if indicator eq 'M' then begin
    
    if (first gt 0) then Frequency = [[[Frequency]],[[dummy]]]
    first=first+1
    dummy=replicate(!values.F_NAN,[64,32])
    
    u=''
    reads,line,y,m,d,h,mm,s,U,dvs,dsn,bw1,cc,mdq,$
      format='(4x,I2,I2,I2,I2,I2,I2,1x,A3,4x,F6,4x,I10,3x,I7,3x,I7,4x,I6)';,4x,I6,4x,E6,4x,A8,4x,F5,4x,I11,3x,I8,4x,I4)'    
    if (first eq 0) then begin
      year=y & month=m & day=d & hour=h & minute=mm & seconds=s & UTC=U
      ;avg=ave & resolution=stp & asl=seal & sample_rate=smp
      ;Noise_0=nf0 & Noise_1=nf1
      ;Version_service=svs
      Version_firmware=dvs
      serial=dsn & bw=bw1
      calibration = cc & valid_spectra=mdq          
    endif else begin
      year=[year,y] & month=[month,m] & day=[day,d] & hour=[hour,h] & minute=[minute,mm] & seconds=[seconds,s] & UTC=[UTC,U]
      ;avg=[avg,ave] & resolution=[resolution,stp] & asl=[asl,seal] & sample_rate=[sample_rate,smp]
      ;Noise_0=[Noise_0,nf0] & Noise_1=[Noise_1,nf1]
      ;Version_service=[Version_service,svs] & 
      Version_firmware=[Version_firmware,dvs]
      serial=[serial,dsn] & bw=[bw,bw1]
      calibration=[calibration,cc] & valid_spectra=[valid_spectra,mdq] 
    endelse  
    ;print,y,m,d,h,mm,s,U,dvs,dsn,bw,cc,mdq  
  ;*** H ***MEASURING HEIGHT  s
  endif else if indicator eq 'H' then begin
    he=fltarr(32)
    reads,line,he,format='(3x,32(I9))'
    if (first eq 0) then begin
      height = he
    endif else begin
      height = [[height],[he]]
    endelse  
  ;*** T ***TRANSFER FUNCTION       
  endif else if indicator eq 'T' then begin
    trans=fltarr(32)
    reads,line,trans,format='(3x,32(F9))'
    if (first eq 0) then begin
      Transfer_function = trans
    endif else begin
      Transfer_function = [[Transfer_function],[trans]]       
    endelse 
  ;*** F ***   64 lines of SPECTRAL FREQUENCY          
  endif else if indicator eq 'F' then begin
    u=0
    reads,line,u,format='(1x,I2)'
    ;print,u
    num=strlen(line)
    num=fix((num-3.)/9.)
    if (num gt 0) then freq=replicate(!values.f_nan,num) else freq=!values.f_nan
    if (num gt 0) then reads,line,freq,format='(3x,32(F9))'
    nan_replace=where(freq eq 0.0,count)
    if (count gt 0) then freq[nan_replace]=!values.F_NAN
    ;if (num le 30) then freq=[freq,nans[num:30]]
    ;if (num eq 0) then freq=nans
    if (first eq 0) then begin
      Frequency[u,*]=freq
    endif else begin
      dummy[u,*]=freq       
    endelse  
  endif    
;endwhile
endfor
free_lun,lun

rot_num=6
julian_day=julday(month,day,year,hour,minute,seconds)
height=rotate(height,rot_num)

;*********************************************************
;  Plot data
;*********************************************************

;  Color common
common colors,r_orig,g_orig,b_orig,r_curr,g_curr,b_curr

;  Set up the device
set_plot,'z'

;  xsize of the plot window
xwin=800  

;  ysize of the plot window
ywin=900  

;  Size the window
device,set_resolution=[xwin,ywin]

;  Load the color table
loadct,39   

;  Make the background white
!p.background=!d.n_colors-1  

top_color=254
;  Date label format
dummy=label_date(date_format=['%H:%I'])

;  Xticks point out from plot
!x.ticklen=-0.02
!y.ticklen=-0.02

;  String positioning
;dr=0.09  ;string offset from the right xcoor
;dl=0.025  ;string offset from the left xcoor
;dt=0.025  ;string offset from the top ycoor

;***********************************************************
;  Plot spectra
;***************************************************
erase
xl=0.1  &  xr=0.85
yb=0.1  &  yt=0.8
sx=0.20
sy=0.15
numplots_x=2  ;number of plots in x direction
numplots_y=2  ;number of plots in y direction

position_plots,xl,xr,yb,yt,sx,sy,numplots_x,numplots_y,pos
cbpos=pos
cbpos[*,0]=pos[*,2]+0.10
cbpos[*,2]=cbpos[*,0]+0.01

dmax=50;max(frequency)
dmin=min(frequency)

dopvel=indgen(64)*0.1905

for i=0,3 do begin

  freqt=frequency[*,*,i]

  pnum=i
  image_var=bytscl(freqt,min=min,max=dmax,top=top_color)

  ;  Set up the plot space
  s=size(freqt,/dimensions)
  contour,freqt,indgen(s[0]),indgen(s[1]),/nodata,$
    xstyle=4,ystyle=4,position=pos[pnum,*],/noerase

  ;  Grab the plot location
  px=!x.window*!d.x_vsize
  py=!y.window*!d.y_vsize
  sx=px[1]-px[0]+1
  sy=py[1]-py[0]+1

  ;  Put up the image
  tv,congrid(image_var,sx,sy),px[0],py[0]

  ;  Put up the axes
  contour,freqt,dopvel,height[0,*]/1000.0,xstyle=1,ystyle=1,/nodata,$
    color=0,$
    ;xtickunits='time',xtickformat='label_date',$
    xtitle='doppler velocity m/s',$;spectral bin',$
    ytitle='height km',$;'range gate',$
    title='20'+string(year[i],format='(I02)')+string(month[i],format='(I02)')+string(day[i],format='(I02)')$
    +' '+string(hour[i],format='(I02)')+':'+string(minute[i],format='(I02)')+':'+string(seconds[i],format='(I02)'),$
    position=pos[pnum,*],/noerase
  ;xyouts,pos[pnum,0]+dl,pos[pnum,3]-dt,'fall_velocity',$
  ;  color=255,/normal
  ;  Color bar
  colorbar_fanning,maxrange=dmax,minrange=dmin,$
    ncolors=top_color,format='(F8.2)',$
    vertical=1,color=0,position=cbpos[pnum,*]
endfor
write_gif,'mrr.raw.gif',tvrd()

stop

end