;IDL program to create plots of along-track integrated Poynting flux for mp/qt at polar cap(pc)/auroral zone/total.
;Sx0 poynting flux data: IDM
;(1) exclude idmqual 3 and 4
;(2) exclude NaN
;Overplot of Poynting Flux


;Define bin size: 2.5 MLat*1MLT
mlat_window = 2.0
binsz = 180./mlat_window ;2.0 MLat
mlat_bin = (findgen(binsz)-(binsz/2.-1.))*mlat_window
mlt_bin = findgen(24)

binave_poynt  = dblarr(24, binsz);*!Values.F_NaN
bin_poynt = dblarr(24,binsz,10000L)*!Values.F_NaN
nbin_poynt = dblarr(24,binsz)
Vsc = 7415. ; m/sec 
min_cb = 50.
max_cb = 100.

;read in DMSP Poynting flux
dir = './moderate_st/';'./strong_st/';'./moderate_st/';'./weak_st/';
filename = file_search(dir + 'Poynt*_mp_n.sav') ;(dir + 'poynt_F16-2011218002250.sav')     
nfiles = n_elements(filename)
print, 'number of files = ',nfiles

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;mp;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
cxpos_CRB1 = dblarr(nfiles*10)*!Values.F_NaN
cypos_CRB1 = dblarr(nfiles*10)*!Values.F_NaN
cxpos_CRB2 = dblarr(nfiles*10)*!Values.F_NaN
cypos_CRB2 = dblarr(nfiles*10)*!Values.F_NaN

ctotal_poyntn = fltarr(3000000L)*!Values.F_NaN
ctotal_mlatn = fltarr(3000000L)*!Values.F_NaN
ctotal_mltn = fltarr(3000000L)*!Values.F_NaN
nptsn0 = 0

for jj = 0, nfiles-1 do begin ;begin major analysis loop, 1 file/loop.
    restore, filename(jj)
    
    indexn = where(finite(total_poyntn),nptsn) ;mp
    
    if (nptsn gt 0) then begin
 
        ctotal_poyntn(nptsn0:nptsn0+nptsn-1) = total_poyntn(indexn)
        ctotal_mlatn(nptsn0:nptsn0+nptsn-1) = total_mlatn(indexn)
        ctotal_mltn(nptsn0:nptsn0+nptsn-1) = total_mltn(indexn)
        print, 'nptsn0=',nptsn0,'    nptsn=',nptsn,'    nptsn0+nptsn=',nptsn0+nptsn
        nptsn0 = nptsn0+nptsn

    endif
endfor                          ;end major analysis loop.


;---------------------Figure(mp)----------------------------------------------------------------
set_plot,'ps'
;device,/inches,xsize=10.90,ysize=5.347,yoffset=10.9,xoffset=1.0,/color,/landscape
device,filename='Poynt' + 'overplot_north_mp.ps',/portrait
device,/color

loadct, 39
; Plot Northern
polar_dial_n,1,1
xyouts,-49.0,44.0,'Poynting Flux(mW/m!U2!N)',charsize=1.1
xyouts,-49.0,40.0,'Northern Hemisphere',charsize=1.1

;sort the poyntn
fp_sort = sort(ctotal_poyntn)

        radius = (90.-ctotal_mlatn(fp_sort))
        angle  = ((ctotal_mltn(fp_sort)-6.)*!pi/12.)
        xpos = radius*cos(angle)
        ypos = radius*sin(angle)
        
 ;;;poynting flux color contour
        loadct, 33
        dcolor = (ctotal_poyntn(fp_sort)-min_cb)/(max_cb-min_cb)*256. ; [0, 50]] 
        for i = 0L, n_elements(ctotal_poyntn)-1 do begin
            usersym,[-1,1,1,-1,-1],[1,1,-1,-1,1],/fill,color = dcolor(i) ;filled square
            ;oplot,xpos(i:i+1),ypos(i:i+1),psym = 8, symsize = 0.8;,thick =20
            plots,xpos(i),ypos(i),psym = 8, symsize = 0.25;0.5
        endfor

;;CRBs
;loadct, 10
;for jj = 0, nfiles-1 do begin ;begin major analysis loop, 1 file/loop.
;    restore, filename(jj)
;    
;    fp_nan  = where(finite(xpos_CRB1), count1)
;    usersym,[-2./sqrt(3),0,2./sqrt(3),-2./sqrt(3)],[-1,1,-1,-1],/fill, color = 200 ;pink filled triangle
;    for inan = 0, count1-1 do begin
;        plots, xpos_CRB1(fp_nan(inan)), ypos_CRB1(fp_nan(inan)), psym = 8, symsize = .5
;    endfor
;    fp_nan  = where(finite(xpos_CRB2), count2)
;    usersym,[-1,1,1,-1,-1],[1,1,-1,-1,1],/fill,color = 200 ;pink filles square
;    for inan = 0, count2-1 do begin
;        plots, xpos_CRB2(fp_nan(inan)), ypos_CRB2(fp_nan(inan)), psym = 8, symsize = .5
;    endfor
;endfor

loadct, 33
;;colorbar
    number_levels = 33
    levels_contour = dblarr(number_levels)
    colors_contour = dblarr(number_levels)
    color_table_offset =.0
maxData = max_cb; max(alog10(temps),/nan)
minData = min_cb                    ; min(alog10(temps),/nan)
for i = 0, number_levels-1 do begin
    levels_contour(i) = (maxData-minData)*i/(number_levels-1.0)+minData
    colors_contour(i) = (255.0-color_table_offset)*i/(number_levels-1.0)+color_table_offset
endfor
contour, [[levels_contour],[levels_contour]],levels_contour,indgen(2),charsize=0.8,charthick=1, $
  levels=levels_contour,c_color=colors_contour,/fill,position=[0.55,0.05,0.85,0.08],/noerase,ystyle=4, $
  xr=[minData,maxData],xstyle=1,xticks = 5

device,/close
set_plot,'x'

;stop

;;;;;;;;;binplots;;;;;;;;;;;;;
fp_ctotal = where(finite(ctotal_poyntn),nptsn) 
ctotal_poyntn = ctotal_poyntn(fp_ctotal)
ctotal_mlatn = ctotal_mlatn(fp_ctotal)
ctotal_mltn = ctotal_mltn(fp_ctotal)

for ipoint = 0L, nptsn-1 do begin
    index_mlat = max(where(ctotal_mlatn(ipoint) ge mlat_bin))
    index_mlt = max(where(ctotal_mltn(ipoint) ge mlt_bin))
    
    bin_poynt(index_mlt,index_mlat,nbin_poynt(index_mlt,index_mlat)) = ctotal_poyntn(ipoint)
    nbin_poynt(index_mlt,index_mlat) = nbin_poynt(index_mlt,index_mlat)+1
    ;print, nbin_poynt(index_mlt,index_mlat)
endfor

set_plot,'ps'
device,filename='Poynt' + strmid(dir,2,11) + '_north_avebin_mp_bin2.ps',/portrait
device,/color
loadct, 39
polar_dial_n,1,1
;xyouts,.0,54.0,title1,alignment=0.5,charsize=1.2
xyouts,-49.0,44.0,'Poynting Flux(mW/m!U2!N)',charsize=1.1
xyouts,-49.0,40.0,'Northern Hemisphere',charsize=1.1

loadct, 33    ;rainbow 0=black, 255=white, 220ish=red, 50ish=blue
binave_poynt = mean(bin_poynt,dimension = 3,/nan) ;max(bin_poynt,dimension =3,/nan);
bin_counts = nbin_poynt
min_binave_poynt = min(binave_poynt,/nan)
max_binave_poynt = max(binave_poynt,/nan)
xyouts,-49.0,-44.0,'min:'+string(min_binave_poynt, format = '(f6.2)'),charsize=1.1
xyouts,-49.0,-48.0,'max:'+string(max_binave_poynt, format = '(f6.2)'),charsize=1.1

;;;plot binned average
bin_color = binave_poynt/20.*256      ;(max(binave_poynt)-min(binave_poynt))*256 ;poynting flux

for imlat = 0, (binsz-2) do begin
    if (mlat_bin(imlat) ge 50.) then begin
        for imlt = 0, 22 do begin
            bin_radius = (90.-mlat_bin(imlat))
            bin_angle  = (mlt_bin(imlt)-6.)*!pi/12.
            bin_xpos1 = bin_radius*cos(bin_angle)
            bin_ypos1 = bin_radius*sin(bin_angle)
            bin_radius = (90.-mlat_bin(imlat+1))
            bin_angle  = (mlt_bin(imlt)-6.)*!pi/12.
            bin_xpos2 = bin_radius*cos(bin_angle)
            bin_ypos2 = bin_radius*sin(bin_angle)
            bin_radius = (90.-mlat_bin(imlat+1))
            bin_angle  = (mlt_bin(imlt+1)-6.)*!pi/12.
            bin_xpos3 = bin_radius*cos(bin_angle)
            bin_ypos3 = bin_radius*sin(bin_angle)
            bin_radius = (90.-mlat_bin(imlat))
            bin_angle  = (mlt_bin(imlt+1)-6.)*!pi/12.
            bin_xpos4 = bin_radius*cos(bin_angle)
            bin_ypos4 = bin_radius*sin(bin_angle)
            xx = [bin_xpos1,bin_xpos2,bin_xpos3,bin_xpos4]
            yy = [bin_ypos1,bin_ypos2,bin_ypos3,bin_ypos4]
            
            if ((bin_counts(imlt,imlat) ne 0.) and finite(binave_poynt(imlt,imlat))) then begin
                polyfill,xx,yy,color = bin_color(imlt,imlat);,/line_fill
            endif else begin
                ;;;missing grids = white
                loadct, 39
                polyfill,xx,yy,color = 255;white
                loadct, 33
            endelse
        endfor
;last mlt_bin bin: 23:24mlt

        bin_radius = (90.-mlat_bin(imlat))
        bin_angle  = (23.-6.)*!pi/12.
        bin_xpos1 = bin_radius*cos(bin_angle)
        bin_ypos1 = bin_radius*sin(bin_angle)
        bin_radius = (90.-mlat_bin(imlat+1))
        bin_angle  = (23.-6.)*!pi/12.
        bin_xpos2 = bin_radius*cos(bin_angle)
        bin_ypos2 = bin_radius*sin(bin_angle)
        bin_radius = (90.-mlat_bin(imlat+1))
        bin_angle  = (24.-6.)*!pi/12.
        bin_xpos3 = bin_radius*cos(bin_angle)
        bin_ypos3 = bin_radius*sin(bin_angle)
        bin_radius = (90.-mlat_bin(imlat))
        bin_angle  = (24.-6.)*!pi/12.
        bin_xpos4 = bin_radius*cos(bin_angle)
        bin_ypos4 = bin_radius*sin(bin_angle)
        xx = [bin_xpos1,bin_xpos2,bin_xpos3,bin_xpos4]
        yy = [bin_ypos1,bin_ypos2,bin_ypos3,bin_ypos4]
        
        if ((bin_counts(23,imlat) ne 0.) and finite(binave_poynt(23,imlat))) then begin
            polyfill,xx,yy,color = bin_color(23,imlat)
        endif else begin
                ;;;missing grids = white
            loadct, 39
            polyfill,xx,yy,color = 255 ;white
            loadct, 33
        endelse
    endif    
endfor

loadct, 39
;additional parts
pi = 3.141592654
angle = 2.*pi*indgen(721)/720.
for i=10.,40.,10. do begin
    radius = replicate(i,721)
    oplot,/polar,radius,angle,linestyle =2
endfor
oplot,[0.,0.],[-45.,45.]
oplot,[-45.,45.],[0.,0.]
for i=10,40,10 do begin
    tang = -45.*pi/180.
    
    txt = ' '+strtrim(string(90-i),2) + '!9%!5' ;northern hemisphere
;   txt = ' '+strtrim(string(i-90),2) + '!9%!5'
    
    if (i eq 40) then txt = txt + ' MLAT'
    xyouts,float(i)*cos(tang),float(i)*sin(tang),txt,charsize=csize
endfor
xyouts,0.,-47.5,'0!CMLT',charsize=csize,alignment=0.5
xyouts,46.,0.,'6',charsize=csize,alignment=0.
xyouts,0.,46.,'12',charsize=csize,alignment=0.5
xyouts,-46.,0.,'18',charsize=csize,alignment=1.

;;;sampled bin: mlt=[12, 13] and mlat=[82,84]
;;sampled bin: mlt=[10, 11] and mlat=[78,80]
loadct, 10
sample_angle = ((10.-6.)+indgen(11)/10.)*!pi/12.
oplot,/polar,replicate((90-80),11),sample_angle,color = 200,thick = 8;84MLAT arc
oplot,/polar,replicate((90-78),11),sample_angle,color = 200,thick = 8;82MLAT arc
sample_radius = (90.-80.)+indgen(21)/10.
oplot,/polar,sample_radius, replicate((10-6),21)*!pi/12.,color = 200,thick = 8;12MLT line
oplot,/polar,sample_radius, replicate((11-6),21)*!pi/12.,color = 200,thick = 8;13MLT line


loadct, 33
;;colorbar
number_levels = 33
levels_contour = dblarr(number_levels)
colors_contour = dblarr(number_levels)
color_table_offset =.0
maxData = 20. ;200.;max(binave_poynt); max(alog10(temps),/nan)
minData = 0.; min(alog10(temps),/nan)
for i = 0, number_levels-1 do begin
   levels_contour(i) = (maxData-minData)*i/(number_levels-1.0)+minData
   colors_contour(i) = (255.0-color_table_offset)*i/(number_levels-1.0)+color_table_offset
endfor
contour, [[levels_contour],[levels_contour]],levels_contour,indgen(2),charsize=0.8,charthick=1, $
levels=levels_contour,c_color=colors_contour,/fill,position=[0.55,0.05,0.85,0.08],/noerase,ystyle=4, $
xr=[minData,maxData],xstyle=1,xticks = 5

loadct, 39
;;;;;;;histogram plot
!p.multi = [0,3,4]
for iimlt = 0, 23 do begin
    for iimlat = (binsz-12), (binsz-1) do begin
        if (total(finite(bin_poynt(iimlt,iimlat,*))) ne 0.) then begin
            hist_mp = histogram(bin_poynt(iimlt,iimlat,*),binsize = 1,locations = binvals,/nan)
            plot,binvals,hist_mp,xtitle = 'Sx(mW/m!U2!N)',ytitle = 'counts',xr=[-20,100],yr=[0,100],psym=-4,symsize = 0.5
            ave_sx = mean(bin_poynt(iimlt,iimlat,*),/nan)
            fp_over = where(binvals ge ave_sx)
            npt_over = total(hist_mp(fp_over))
            npt_total = total(hist_mp)
            ;print, 'over=',npt_over,'   total=',npt_total
            oplot, [ave_sx,ave_sx],[0, 100],color= 250
            
            xyouts, 10,60, 'mlt='+string(mlt_bin(iimlt),format='(f4.1)')+ ' mlat='+string(mlat_bin(iimlat),format='(f5.1)')+'!9%!5',charsize = 0.5
            xyouts, 10,50, 'average='+string(mean(bin_poynt(iimlt,iimlat,*),/nan),format='(f5.2)')+' mW/m!U2!N',charsize=0.5
            xyouts, 10, 40, 'over#/total#='+string(npt_over,format='(i5)')+'/'+string(npt_total,format='(i5)'),charsize=0.5
;hist_mp = histogram(bin_poynt(iimlt,iimlat,*),binsize = 5,locations = binvals,/nan)
;plot,binvals,hist_mp,xtitle = 'Sx(mW/m!U2!N)',ytitle = 'counts',title =
;'binsize = 5'
        endif else begin
            plot,[-20,40],[1,20],/nodata,xr = [-20,100],yr = [0,100]
            xyouts, 10,50, 'mlt='+string(mlt_bin(iimlt),format='(f4.1)')+ ' mlat='+string(mlat_bin(iimlat),format='(f5.1)'),charsize = 0.5
        endelse
    endfor
endfor
!p.multi = 0

device,/close
set_plot,'x'


;;;;;;save mp Sx (not sorted) data to /hao/aim3/huangys/tiegcm1.94_data/temp/DMSP_PoyntFlux
;save,total_poyntn,total_mlatn,total_mltn,xpos_CRB1,ypos_CRB1,xpos_CRB2,ypos_CRB2, $
;  filename='/hao/aim3/huangys/tiegcm1.94_data/temp/DMSP_PoyntFlux/Poynt'+strmid(dir,2,11)+'_mp_n.sav'

end