; DS3.PRO  -  make figure 3

; set the color table 
  
  device,true=24
  device,de=0
  device,retain=2
;  loadct,39  ; the one mei-ching fok used
;  loadct,41 ;our usual

  nz  = 304
  nf  = 89
  nt  =  200

  dir  = './' 

  day1 = 279  

  close,1
  time = fltarr(5,nt)
  openr,1,dir+'ds2_time.dat'
  print,' reading time'
  readf,1,time
  close,1

  time2 = fltarr(5,nt)
  openr,1,dir+'ds2_time_2s.dat'
  print,' reading time 2s'
  readf,1,time2
  close,1

  day = fltarr(nt)
  day(0) = 279.
  for i = 1,nt-1 do begin
    day(i) = day(i-1)
    if (time(1,i) lt time(1,i-1)) then day(i) = day(i) + 1.
  endfor
  day2 = fltarr(nt)
  day2(0) = 279.
  for i = 1,nt-1 do begin
    day2(i) = day2(i-1)
    if (time2(1,i) lt time2(1,i-1)) then day2(i) = day2(i) + 1.
  endfor

; latitude

  close,1
  glat = fltarr(nz,nf)
  openr,1,dir+'ds2_glat.dat',/f77_unformatted
  print,' reading glat'
  readu,1,glat

; altitude

  close,1
  zalt = fltarr(nz,nf)
  openr,1,dir+'ds2_zalt.dat',/f77_unformatted
  print,' reading zalt'
  readu,1,zalt

; density data

  deni1   = fltarr(nz,nf,nt)
  deni5   = fltarr(nz,nf,nt)
  deni1s = fltarr(nz,nf,nt)
  deni1n = fltarr(nz,nf,nt)
  deni52  = fltarr(nz,nf,nt)
  close,1
  openr,1,dir+'ds3_denh.dat',/f77_unformatted
  print,' reading denh'
  readu,1,deni1
  close,1
  openr,1,dir+'ds3_denhe.dat',/f77_unformatted
  print,' reading denhe'
  readu,1,deni5
  close,1
  openr,1,dir+'ds3_denhs_2s.dat',/f77_unformatted
  print,' reading denhs 2s'
  readu,1,deni1s
  close,1
  openr,1,dir+'ds3_denhn_2s.dat',/f77_unformatted
  print,' reading denhn 2s'
  readu,1,deni1n
  close,1
  openr,1,dir+'ds3_denhe_2s.dat',/f77_unformatted
  print,' reading denhe 2s'
  readu,1,deni52

; velocity data

  vsi1   = fltarr(nz,nf,nt)
  vsi5   = fltarr(nz,nf,nt)
  vsi1s = fltarr(nz,nf,nt)
  vsi1n = fltarr(nz,nf,nt)
  vsi52  = fltarr(nz,nf,nt)
  close,1
  openr,1,dir+'ds3_vh.dat',/f77_unformatted
  print,' reading vh'
  readu,1,vsi1
  close,1
  openr,1,dir+'ds3_vhe.dat',/f77_unformatted
  print,' reading vhe'
  readu,1,vsi5
  close,1
  openr,1,dir+'ds3_vhs_2s.dat',/f77_unformatted
  print,' reading vhs 2s'
  readu,1,vsi1s
  close,1
  openr,1,dir+'ds3_vhn_2s.dat',/f77_unformatted
  print,' reading vhn 2s'
  readu,1,vsi1n
  close,1
  openr,1,dir+'ds3_vhe_2s.dat',/f77_unformatted
  print,' reading vhe 2s'
  readu,1,vsi52

; specity indices and arrays to be plotted
  
  nfm = 70
  ntm = 94
  fna1  = reform(deni1 (*,nfm,*))
  fna2  = reform(deni5 (*,nfm,*))
  fnb2 = reform(deni52(*,nfm,*))
  fnb3 = reform(deni1s(*,nfm,*))
  fnb4 = reform(deni1n(*,nfm,*))  
  fnb1 = fnb3 + fnb4
  yt = '!N[H!U+!N], [He!U+!N] (cm!U-3!N)'
  yr = [0.05,30.]
  ytp = 1
  fnc1  = reform(vsi1 (*,nfm,*))/1.e5
  fnc2  = reform(vsi5 (*,nfm,*))/1.e5
  fnd2 = reform(vsi52(*,nfm,*))/1.e5
  fnd3 = reform(vsi1s(*,nfm,*))/1.e5
  fnd4 = reform(vsi1n(*,nfm,*))/1.e5
; density-weighted average  
  fnd1 = (fnd3*fnb3 + fnd4*fnb4)/fnb1
  yt2 = '!Nv!BH!N+, v!BHe!N+ (km/s)'
  yr2 = [-10,10]
  ytp2 = 0


   loadct,10   ;green blue pink

;  loadct,41 ;our usual

  re = 6370. 
; -----------------------
; set plot to the screen
; -----------------------

  set_plot,'x'

; set to a single window

  !p.multi=0

; set the background to white 
; (instead of the default which is black)

  !p.background=65535

; set the window size

  xwin = 600 
  ywin = 600
  window,xsize=xwin,ysize=ywin

; scale the window for the postscript plot

  xps  = 6.0
  yps  = xps * ywin / xwin

; positions
  pos1 = [.16,.53,.56,.92]
  pos3 = [.16,.12,.56,.51]
  pos2 = [.59,.53,.99,.92]
  pos4 = [.59,.12,.99,.51]

; position and size of the color bar

  posc = [.40,.98,.70,.98]
  xszc = ( posc(2) - posc(0) ) * xps
  yszc = ( posc(3) - posc(1) ) * yps
  xstc = posc(0) * xps
  ystc = posc(1) * yps
  xszw = ( posc(2) - posc(0) ) * xwin

; set font to simplex roman

  xyouts,.5,.5,'!3 ',/normal 

  xt = '!NGeographic Latitude'
  samis  = '!NSAMI3 (One fluid H!U+!N)'
  samiss = '!NSAMI3 (Two fluid H!U+!N)'
  
; set specific minimum and maximun values

;  fmn    =  0.
;  fmx    =  1.e4
  fmn    = -2.
  fmx    = 5.
  glatmn  = -64.
  glatmx  = 62.
  zaltmn  = 0.  ; in Re
  zaltmx  = 4.5
  
    hr   = time(1,ntm)
    min  = time(2,ntm)
    if min ge 59.0 then hr = hr+1
    if min ge 59.0 then min = 0.
    phr  = string(format='(i2)',hr)
    while(((i=strpos(phr,' '))) ne -1) do strput,phr,'0',i
    pmin = string(format='(i2)',min)
    while(((i=strpos(pmin,' '))) ne -1) do strput,pmin,'0',i
    days = 'Day '+string(format='(i3)',day(ntm))
    uts  = phr+pmin+' UT' 

    hr   = time2(1,ntm)
    min  = time2(2,ntm)
    if min ge 59.0 then hr = hr+1
    if min ge 59.0 then min = 0.
    phr  = string(format='(i2)',hr)
    while(((i=strpos(phr,' '))) ne -1) do strput,phr,'0',i
    pmin = string(format='(i2)',min)
    while(((i=strpos(pmin,' '))) ne -1) do strput,pmin,'0',i
    days2 = 'Day '+string(format='(i3)',day2(ntm))
    uts2    = days2+'  '+phr+pmin+' UT';+days


  print,' uts  = ',uts
  print,' uts2 = ',uts2
  
; label the L value of the field line
  al = zalt(nz/2,nfm)/re + 1.0
  als = 'L = '+string(format='(f3.1)',al)
  print,als
  
; plot the lines
 
  lst2 = 2
  lcol2 = 0
  ics = 130  ; 130    60 (ct 41)
  icn =  25  ;  25   150 (ct 41)
  plot,glat(*,nfm),fna1(*,ntm),$
          yrange=yr,ystyle=1,ytype=ytp, $
          xrange=[glatmn,glatmx],xstyle=1,$
          xticknam=[' ',' ',' ',' ',' ',' ',' ',' '], $
          xtitle=' ',$
          ytitle=yt,$
          color=0,pos=pos1
  oplot,glat(*,nfm),fna2(*,ntm),linestyle=lst2,color=lcol2

  plot,glat(*,nfm),fnb1(*,ntm),$
          yrange=yr,ystyle=1,ytype=ytp, $
          xrange=[glatmn,glatmx],xstyle=1,$
          xticknam=[' ',' ',' ',' ',' ',' ',' ',' '], $
          yticknam=[' ',' ',' ',' ',' ',' ',' ',' '], $
          xtitle=' ',$
          ytitle=' ',$
          color=0,pos=pos2,/noerase
  oplot,glat(*,nfm),fnb2(*,ntm),linestyle=lst2,color=lcol2
  oplot,glat(*,nfm),fnb3(*,ntm),linestyle=1,color=ics
  oplot,glat(*,nfm),fnb4(*,ntm),linestyle=1,color=icn
  
  plot,glat(*,nfm),fnc1(*,ntm), $
          yrange=yr2,ystyle=1,ytype=ytp2, $
          xrange=[glatmn,glatmx],xstyle=1,$
 ;         xticknam=[' ',' ',' ',' ',' ',' ',' ',' '], $
          xtitle=xt,$
          ytitle=yt2,$
          color=0,pos=pos3,/noerase
  oplot,glat(*,nfm),fnc2(*,ntm),linestyle=lst2,color=lcol2

  plot,glat(*,nfm),fnd1(*,ntm), $
          yrange=yr2,ystyle=1,ytype=ytp2, $
          xrange=[glatmn,glatmx],xstyle=1,$
;          xticknam=[' ',' ',' ',' ',' ',' ',' ',' '], $
          yticknam=[' ',' ',' ',' ',' ',' ',' ',' '], $
          xtitle=xt,$
          ytitle=' ',$
          color=0,pos=pos4,/noerase
  oplot,glat(*,nfm),fnd2(*,ntm),linestyle=lst2,color=lcol2
  oplot,glat(*,nfm),fnd3(*,ntm),linestyle=1,color=ics
  oplot,glat(*,nfm),fnd4(*,ntm),linestyle=1,color=icn
  
; print the times

  xp1 = .18
  xp2 = .61
  xp3 = .50
  xp4 = .93
  yp1 = .265
  yp2 = .475
  yp3 = .685
  yp4 = .885
  yp5 = .935

  sz = 1.2
  xyouts,xp2-.015,yp5,uts,size=sz,color=0,/normal
  xyouts,xp2-.015,yp5+0.03,days,size=sz,color=0,/normal
  
  xyouts,xp1-0.01,yp5,samis,color=0,/normal
  xyouts,xp2+0.13,yp5,samiss,color=0,/normal
  xyouts,xp3-0.04,yp5,als,color=0,/normal
  xyouts,0.70,0.70,'!3S',size=sz,color=ics,/normal
  xyouts,0.85,0.70,'!3N',size=sz,color=icn,/normal
  xyouts,0.72,0.47,'!3S',size=sz,color=ics,/normal
  xyouts,0.83,0.17,'!3N',size=sz,color=icn,/normal
  
  xyouts,xp3,yp4,'!3(a)',size=sz,color=0,/normal
  xyouts,xp4,yp4,'!3(c)',size=sz,color=0,/normal
  xyouts,xp3,yp2,'!3(b)',size=sz,color=0,/normal
  xyouts,xp4,yp2,'!3(d)',size=sz,color=0,/normal
  
; -----------------------
; set plot to postscript
; -----------------------

  !p.thick=4
  !p.charthick=4
  !x.thick=4
  !y.thick=4
  !p.charsize = 1.0

  set_plot,'ps'
  device,file='fig.ps',bits_per_pixel=8,/color, $
         xsize=xps,ysize=yps,/inches,xoffset=1.25,yoffset=1.25

; plot the line
  plot,glat(*,nfm),fna1(*,ntm),$
          yrange=yr,ystyle=1,ytype=ytp, $
          xrange=[glatmn,glatmx],xstyle=1,$
          xticknam=[' ',' ',' ',' ',' ',' ',' ',' '], $
          xtitle=' ',$
          ytitle=yt,$
          color=0,pos=pos1
  oplot,glat(*,nfm),fna2(*,ntm),linestyle=lst2,color=lcol2

  plot,glat(*,nfm),fnb1(*,ntm),$
          yrange=yr,ystyle=1,ytype=ytp, $
          xrange=[glatmn,glatmx],xstyle=1,$
          xticknam=[' ',' ',' ',' ',' ',' ',' ',' '], $
          yticknam=[' ',' ',' ',' ',' ',' ',' ',' '], $
          xtitle=' ',$
          ytitle=' ',$
          color=0,pos=pos2,/noerase
  oplot,glat(*,nfm),fnb2(*,ntm),linestyle=lst2,color=lcol2
  oplot,glat(*,nfm),fnb3(*,ntm),linestyle=1,color=ics
  oplot,glat(*,nfm),fnb4(*,ntm),linestyle=1,color=icn

  plot,glat(*,nfm),fnc1(*,ntm), $
          yrange=yr2,ystyle=1,ytype=ytp2, $
          xrange=[glatmn,glatmx],xstyle=1,$
;          xticknam=[' ',' ',' ',' ',' ',' ',' ',' '], $
          xtitle=xt,$
          ytitle=yt2,$
          color=0,pos=pos3,/noerase
  oplot,glat(*,nfm),fnc2(*,ntm),linestyle=lst2,color=lcol2
  
  plot,glat(*,nfm),fnd1(*,ntm), $
          yrange=yr2,ystyle=1,ytype=ytp2, $
          xrange=[glatmn,glatmx],xstyle=1,$
;          xticknam=[' ',' ',' ',' ',' ',' ',' ',' '], $
          yticknam=[' ',' ',' ',' ',' ',' ',' ',' '], $
          xtitle=xt,$
          ytitle=' ',$
          color=0,pos=pos4,/noerase
  oplot,glat(*,nfm),fnd2(*,ntm),linestyle=lst2,color=lcol2
  oplot,glat(*,nfm),fnd3(*,ntm),linestyle=1,color=ics
  oplot,glat(*,nfm),fnd4(*,ntm),linestyle=1,color=icn
   
; print the labels

  sz = 1.0
  ush = 0.003
  xyouts,xp2-.015,yp5,uts,size=sz,color=0,/normal
  xyouts,xp2-.015,yp5+0.03,days,size=sz,color=0,/normal

  xyouts,xp1-0.01,yp5,samis,color=0,/normal
  xyouts,xp2+0.13,yp5,samiss,color=0,/normal
  xyouts,xp3-0.04,yp5,als,color=0,/normal
  xyouts,0.70,0.70,'!3S',size=sz,color=ics,/normal
  xyouts,0.85,0.70,'!3N',size=sz,color=icn,/normal
  xyouts,0.72,0.47,'!3S',size=sz,color=ics,/normal
  xyouts,0.83,0.17,'!3N',size=sz,color=icn,/normal

  xyouts,xp3,yp4,'!3(a)',size=sz,color=0,/normal
  xyouts,xp4,yp4,'!3(c)',size=sz,color=0,/normal
  xyouts,xp3,yp2,'!3(b)',size=sz,color=0,/normal
  xyouts,xp4,yp2,'!3(d)',size=sz,color=0,/normal

; close the postscript file 'contour.ps'

  device,/close

; reset the display to the screen

  set_plot,'x'

  !p.charthick = 1
  !p.thick     = 1
  !x.thick=1
  !y.thick=1
  !p.charsize = 1.2

end