;IDL program to create polar plots of CHAMP density data mapped
;to 360km using MSIS.

;time interval for the main phase
time_mp1 = 162.0;0.0;162.;21.60;7.3
;0.0;217.6;215;217.6;22.1;269.4;297.6;217.6;69.0479;67.1840
;(68.4625;67.000);22.24;297.800;269.526;217.796  ;unit:day
time_mp2 = 163.0;366.;164.;22.40;8.3
;366.;218.4;223;218.4;23.6;270.2;298.5;218.4;69.3424;67.3861
;(69.3438;67.625);23.22;298.055;269.889;218.144
mp1 = 0.0;21.7243; 7.39375;plot vertical onset time
;make sure to change the year for read in data file!!!!!!

abs_cut1 = 1.e-15;1.5e-15(Oct);1.e-15(Sep);.6e-15(Aug);1.0e-15(Jan2012)
abs_cut2 = 2.5e-15;2.e-14;3.e-14;1.5e-14;2.e-14;polar plot limits
percent_cut = 30.;30.(Oct);35(Sep);40.(Aug);30.(Jan2012)
plot_cut = 0.0;0.0;2.0e-15;1.5e-15;1.0e-15
Vsc = 7415. ; m/sec 

start_end_year = 2009;indgen(10)+2001

for jj = 0, n_elements(start_end_year)-1 do begin ;begin major analysis loop, 1 file/loop.
print, start_end_year(jj)
;;;;;read in omni2 minutely By and Bz----------------------------------------------
omni_minute = read_omni_minutely_ByBz('../omni_min_ByBz_'+strtrim(start_end_year(jj),2)+'.lst')
year_omni = omni_minute.year
doy_omni = omni_minute.doy
hour_omni = omni_minute.hour
minute_omni = omni_minute.minute

time_omni = double(doy_omni) + (hour_omni+minute_omni/60.)/24. ;unit:day

imf = omni_minute.imf
by = omni_minute.by
bz = omni_minute.bz
V_sw = omni_minute.V_sw
Pres_sw = omni_minute.Pres_sw
symh = omni_minute.sym_h

;remove bad data points
fp = where(imf ge 9999.99)
imf(fp) = !Values.F_NAN
fp = where(by ge 9999.99)
by(fp) = !Values.F_NAN
fp = where(bz ge 9999.99)
bz(fp) = !Values.F_NAN
fp = where(V_sw ge 99999.9)
V_sw(fp) = !Values.F_NAN
fp = where(Pres_sw ge 99.99)
Pres_sw(fp) = !Values.F_NAN
;fp = where(abs(symh) ge 9999.)
;symh(fp) = !Values.F_NAN

;stop

;;;;;read in CHAMP density----------------------------------------------
filename = file_search('./original_data/pre_process/CHAMP_MSIS_400_360_320km_'+strtrim(start_end_year(jj),2)+'_162_highres.dat')

    data = read_champ_msis_orig(filename)
    
    time = data.doy
    alt = data.alt
    lon = data.lon
    lat = data.lat
    stl = data.stl
    ;mlon = data.mlon
    ;mlat = data.mlat
    ;mlt = data.mlt
    den_champ = data.den_champ
    den_msis = data.den_msis
    den_champ_400 = data.den_champ_400
    den_champ_360 = data.den_champ_360
    den_champ_320 = data.den_champ_320

;;;;;;;convert geographic coor to AACGM coor;;;;;;;;;;;;;;;;;;;;;
;IDL>.com aacgmidl.pro
;IDL>.com cnvtime.pro
;IDL>aacgmidl
;IDL>load_coef,"./AACGM_coef/aacgm_coeffs2010.asc"
    nsz = n_elements(time)
    mlat_aacgm = fltarr(nsz)*!Values.F_NaN
    mlon_aacgm = fltarr(nsz)*!Values.F_NaN
    mlt_aacgm = fltarr(nsz)*!Values.F_NaN
    
    year = start_end_year(jj)   ;fix(date/1.e3)
    yrsec = time*3600*24.       ;sec past the start of the year
    for ia = 0, nsz-1 do begin
        cnv_aacgm,lat(ia), lon(ia),alt(ia),mlat_out, mlon_out
        mlat_aacgm(ia) = mlat_out
        mlon_aacgm(ia) = mlon_out
        mlt_aacgm(ia) = calc_mlt(year, yrsec(ia), mlon_out)
    endfor
    
mlat = mlat_aacgm
mlon = mlon_aacgm
mlt = mlt_aacgm

;stop
;;;;;;;interpol omni miniutely data to CHAMP data time step(30~90sec)--------------
by_int = interpol(by, time_omni, time)
bz_int = interpol(bz, time_omni, time)
Vsw_int = interpol(V_sw, time_omni, time)
Pres_int = interpol(Pres_sw, time_omni, time)
symh_int = interpol(symh, time_omni, time)
        
;;;;;;;;;;;;;;;;;;;;;;;;;;;;CHAMP density for mp;;;;;;;;;;;;;;;;;;;;;;;;;
    indexmp = where((time ge time_mp1) and (time le time_mp2),nptsmp) ;mp
    if (nptsmp gt 0) then begin
        timemp = time(indexmp)
        glatmp = lat(indexmp)
        glonmp = lon(indexmp)
        gltmp = stl(indexmp)
        mlatmp = mlat(indexmp)
        mlonmp = mlon(indexmp)
        mltmp = mlt(indexmp)
        den_champ_360mp = den_champ_360(indexmp)
        by_intmp = by_int(indexmp)
        bz_intmp = bz_int(indexmp)
        Vsw_intmp = Vsw_int(indexmp)
        Pres_intmp = Pres_int(indexmp)
        symh_intmp = symh_int(indexmp)
        
        for kk = 0L, nptsmp-2 do begin ;Take care of MLT time singularity.
            if (abs(mltmp(kk)-mltmp(kk+1)) gt 12) then begin
                ikk = where(mltmp ge 12)
                mltmp(ikk) = mltmp(ikk)-24.
            endif
        endfor
        ineg = where(mltmp lt 0.,nptneg) ;Make MLT positive
        if(nptneg gt 0) then mltmp(ineg) = mltmp(ineg) + 24.  
    endif

    den_sm = smooth(den_champ_360mp,130,/nan) ;~22min, 90 degee lat
    delta_den = den_champ_360mp - den_sm
    percent = (den_champ_360mp-den_sm)/den_sm*100
;;;;;;;;finding peaks;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    d0 = den_champ_360mp - shift(den_champ_360mp,1);shift to right
    d1 = den_champ_360mp - shift(den_champ_360mp,-1);shift to left
;    fp_peak = where((d0 gt 0) and (d1 gt 0) and (percent ge percent_cut) and (den_champ_360mp ge abs_cut1),npeak)
    fp_peak = where((d0 gt 0) and (d1 gt 0) and (percent ge percent_cut),npeak)

;;;;;;;;;;calculate hemisphric average without peaks(NH/SH)
    den_withpeak = den_champ_360mp
    den_nopeak = den_withpeak
    den_nopeak(fp_peak) = !values.f_nan

    set_plot,'ps'
    device,filename='CHAMP_den_360km_mp.ps',/landscape ;/portrait
    device,/color
    
    loadct, 39
    !p.multi = [0,1,3]
    
    plot, timemp, den_champ_360mp,ystyle = 8,xstyle = 1,ytitle='Den_360km', $
      xr = [time_mp1, time_mp2],charsize = 1.5,xticks = 5,yr=[0,abs_cut2]
    loadct, 18
    ;oplot, [0, 366],[abs_cut1, abs_cut1], color = 30
    loadct, 39
    oplot, [mp1, mp1],[0,1], color = 125, thick = 6
    oplot, timemp,den_sm,color = 250
    oplot, timemp(fp_peak),den_champ_360mp(fp_peak),color = 250, psym = 1, symsize = 1.
    axis,yaxis = 1,yr=[-90,90],ystyle = 1,/save,ytitle = 'MLAT',charsize = 1.5, yticks = 4
    oplot,timemp,mlatmp,color = 50
    oplot,[0, 365],[0,0],linestyle = 2,color = 50
    ;xyouts, timemp(150), 45., 'threshold = '+string(abs_cut1, format = '(e7.1)')

;;;;;delta density
    plot, timemp, delta_den,ystyle = 8,xstyle = 1,ytitle='delta_Den_360km', $
      yr=[-1.e-15,1.e-15], xr = [time_mp1, time_mp2], charsize = 1.5,xticks = 5
    oplot, [mp1, mp1],[-1.,1.], color = 125, thick = 6
    axis,yaxis = 1,yr=[-90,90],ystyle = 1,/save,ytitle = 'MLAT',charsize = 1.5, yticks = 4
    oplot,timemp,mlatmp,color = 50
    oplot,[0, 365],[0,0],linestyle = 2,color = 50
    
;;;;;percentage
    plot, timemp, percent, ystyle = 9,xstyle = 1,ytitle='%Den_360km',xtitle = 'DOY', $
      yr = [-150,150], xr = [time_mp1, time_mp2],charsize = 1.5,xticks = 5
    loadct, 18
    oplot, [0, 366], [percent_cut, percent_cut], color = 30
    loadct, 39
    oplot, [mp1, mp1],[-200,200], color = 125, thick = 6
    oplot, timemp(fp_peak),percent(fp_peak),color = 250, psym = 1, symsize = 1.
    axis,yaxis = 1,yr=[-90,90],ystyle = 1,/save,ytitle = 'MLAT',charsize = 1.5, yticks = 4
    oplot,timemp,mlatmp,color = 50
    oplot,[0, 365],[0,0],linestyle = 2,color = 50
    ;xyouts, timemp(150), 45., 'threshold = ' + strmid(percent_cut,6,2)

    !p.multi = 0
    device,/close
    set_plot,'x'   
endfor                          ;end major analysis loop.

end