pro plot_poynflux

; am 1/04
; basic program to plot poynting flux from DE2 data in polar
; plot
  close,1
;
  datafile='poynt_60.out'  ; file name
  work = '/tmp/'           ; directory
  path = work + datafile   ; path to file
  openr,1,path             ; open file
;
  nmax  = 253658           ; max. no. of lines in file use wc to get it
  nread =  20000           ; lines to be read in
  ndim  = 3                ; poynting flux 
                           ; rpa_poynting_Wm2,idm_poynting_Wm2,tot_poynting_Wm2
  mlon_pol = 71.0
;
  lons = fltarr(nread)     ; longitudes
  lats = fltarr(nread)     ; latitudes
  data = fltarr(nread,ndim); data points
  year = 0
  day  = 0
  pass_time = 0.
  jbin = 0
  lat = 0.
  mlon= 0.
  mlt = 0.
  rpa = 0.
  idm = 0.
  tot = 0.
;	  
  outfile='./plt.ps'       ; file name for plot
  print,'plotfile name: ',outfile
  set_plot, 'ps'
  device, file=outfile,/portrait, $
       /color, font_size=18, /times
  TVLCT, [0,255,0,0], [0,0,255,0], [0,0,0,255]
;
; read in data
  for n=0,nread-1 do begin
    readf,1,format='(i5,i3,4f9.2,1i7,3e12.4)',year,day,pass_time, $
       lat,mlon,mlt,jbin,rpa,idm,tot
;       print,n,year,day,pass_time,lat,mlt,jbin,rpa,idm,tot
;       
       ut = pass_time/3600.	; get ut pass_time [s]
       lats(n) = lat
       lons(n) = mlon           ; get mag. longitude
       data(n,0) = rpa		; rpa_poynting_Wm2
       data(n,1) = idm  	; idm_poynting_Wm2
       data(n,2) = tot  	; tot_poynting_Wm2 
       data(n,2) = sin(lats(n)*3.14/180.)  	; tot_poynting_Wm2 
  endfor
  print,'lat_range: ',min(lats),max(lats)
  print,'lon_range: ',min(lons),max(lons)
;  
  close,1
;
  crt_nh  =  90.           ; center of plot NH
  crt_sh  = -90.           ; center of plot SH
  crt_nh_min  =  40.       ; min of plot NH
  crt_sh_max  = -40.       ; max of plot SH
  crt_lon =   0.           ; center of plot longitude
  mappos  = [0,0,1.,1.]      ; map position
;  
; NH
  map_set,crt_nh,crt_lon,/stereo ,/noborder ,$
;    limit=[crt_nh_min,-180.,crt_nh,180.] ; limit=[Latmin, Lonmin, Latmax, Lonmax]
    limit=[40.,-180.,90.,180.] 
;  
  triangulate,  lons,lats,tri
  contour,data(*,2),lons,lats,triangulation = tri,$
  /noerase,xrange=[-180,180],yrange=[40,90.]
;
device,/close
;
end
;