;IDL program to create polar plots of CHAMP density data mapped ;to 360km using MSIS. ;time interval for the main phase time_mp1 = 0.;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 = 366.;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 = 5.0e-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 = 2006;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('./CHAMP_MSIS_400_360_320km_'+strtrim(start_end_year(jj),2)+'.dat') data = read_champ_msis(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 ;;;;;;;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 ;;;;% of CHAMP orbit reaching polar latitudes >= 75 MLAT dayofyear = findgen(366)+1.5 ; [1.5:1:366.5] per_mlat = dayofyear*!values.f_nan for idoy = 0, n_elements(dayofyear)-1 do begin fp_doy = where((abs(timemp-dayofyear(idoy)) le 0.5), count_doy) fp_mlat = where((abs(timemp-dayofyear(idoy)) le 0.5) and (abs(mlatmp) ge 75.), count_mlat) ;if (count_doy ne 0.) then per_mlat(idoy) = count_mlat/count_doy*100. endfor set_plot, 'ps' device, filename='CHAMP_360km_AACGM_orbit_75MLAT'+strtrim(start_end_year(jj),2)+'.ps',/portrait device,/color plot, dayofyear, per_mlat, xtitle='DOY',ytitle='percent of CHAMP orbit>75MLAT',xstyle=1 device,/close set_plot,'x' den_sm = smooth(den_champ_360mp,30,/nan) 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 ;;seperate NH/SH ; lat_next = timemp*!values.f_nan ; lat_next(0:nptsmp-2) = mlatmp(1:nptsmp-1) ; multi = mlatmp*lat_next ; fp_traj = where(multi le 0.0,ntraj) ; ; orbit_time = fltarr(ntraj+1)*!values.f_nan ; orbit_den_withpeak = fltarr(ntraj+1)*!values.f_nan ;orbit averaged den with peaks ; orbit_den_nopeak = fltarr(ntraj+1)*!values.f_nan ;;for the first orbit ; orbit_time(0) = mean(timemp(0:fp_traj(0)),/nan) ; orbit_den_withpeak(0) = mean(den_withpeak(0:fp_traj(0)),/nan) ; orbit_den_nopeak(0) = mean(den_nopeak(0:fp_traj(0)),/nan) ; for itra = 0L, ntraj-2 do begin ; orbit_time(itra+1) = mean(timemp(fp_traj(itra):fp_traj(itra+1)),/nan) ; orbit_den_withpeak(itra+1) = mean(den_withpeak(fp_traj(itra):fp_traj(itra+1)),/nan) ; orbit_den_nopeak(itra+1) = mean(den_nopeak(fp_traj(itra):fp_traj(itra+1)),/nan) ; endfor ;;for the last orbit ; orbit_time(ntraj) = mean(timemp(fp_traj(ntraj-1):nptsmp-1),/nan) ; orbit_den_withpeak(ntraj) = mean(den_withpeak(fp_traj(ntraj-1):nptsmp-1),/nan) ; orbit_den_nopeak(ntraj) = mean(den_nopeak(fp_traj(ntraj-1):nptsmp-1),/nan) ; ;;;;;orbit average ; norbit = fix((ntraj+1)/2.) ; orbit_time2 = fltarr(norbit)*!values.f_nan ; orbit_den_withpeak2 = fltarr(norbit)*!values.f_nan ; orbit_den_nopeak2 = fltarr(norbit)*!values.f_nan ; for ior = 0L, norbit-1 do begin ; orbit_time2(ior) = (orbit_time(ior*2)+orbit_time(ior*2+1))/2. ; orbit_den_withpeak2(ior) = (orbit_den_withpeak(ior*2)+orbit_den_withpeak(ior*2+1))/2. ; orbit_den_nopeak2(ior) = (orbit_den_nopeak(ior*2)+orbit_den_nopeak(ior*2+1))/2. ; endfor ;;;orbit average (window size = 1 orbit) with 5-minutes moving step. ; npoints = fix((time_mp2-time_mp1)*24.*60./5.)+1 ; OA_time = timemp(0)+94./2/60./24.+findgen(npoints)*5./60./24. ; OA_den = fltarr(npoints)*!values.f_nan ; for ipt = 0, npoints-1 do begin ; fp_orbit = where(abs(timemp-OA_time(ipt)) le (94./2/60./24.), np_orbit) ; if (np_orbit ne 0) then begin ; OA_den(ipt) = mean(den_champ_360mp(fp_orbit),/nan) ; endif ; endfor ; set_plot,'ps' ; device,filename='OA_CHAMP_den_360km_5min_rm.ps',/landscape ;/portrait ; device,/color ; loadct, 39 ; plot, timemp, den_champ_360mp, xr=[216, 222], xtitle='DOY', ytitle='CHAMP_360km(g/cm3)',psym = -3 ; ;oplot, OA_time, OA_den, color = 250, psym = -3 ; ;legend, ['CHAMP','5-min OA CHAMP'],color =[0, 250], psym = [-3,-3] ; device,/close ; set_plot,'x' 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)') ;;;;;orbit averaged density with/without peaks ; plot, orbit_time, orbit_den_withpeak, psym = -1,ystyle = 8,xstyle = 1,ytitle='Orbit_ave_Den_360km', $ ; xr = [time_mp1, time_mp2],charsize = 2.,xticks = 4 ; oplot, orbit_time2,orbit_den_withpeak2,linestyle = 2 ; oplot, orbit_time, orbit_den_nopeak, psym = -2, color = 240 ; oplot, orbit_time2,orbit_den_nopeak2,linestyle = 2,color = 240 ; legend,['HA with peaks','OA with peaks','HA without peaks','OA without peaks'], $ ; psym = [-1,0,-2,0],colors=[0,0,240,240],linestyle=[0,2,0,2] ; axis,yaxis = 1,yr=[-90,90],ystyle = 1,/save,ytitle = 'MLAT',charsize = 2., yticks = 4 ; oplot,timemp,mlatmp,color = 50 ; oplot,[0, 365],[0,0],linestyle = 2,color = 50 ;;;;;;difference between orbit averaged density with/without peaks ; plot, orbit_time, orbit_den_withpeak-orbit_den_nopeak, psym = -1,ystyle = 8,xstyle = 1, $ ; ytitle='diff_Orbit_ave_Den',xr = [time_mp1, time_mp2],charsize = 2.,xticks = 4 ; oplot, orbit_time2, orbit_den_withpeak2-orbit_den_nopeak2,linestyle=2, color = 240 ; legend,['HA','OA'],psym = [-1,0],color=[0, 200],linestyle=[0,2] ; axis,yaxis = 1,yr=[-90,90],ystyle = 1,/save,ytitle = 'MLAT',charsize = 2., yticks = 4 ; oplot,timemp,mlatmp,color = 50 ; oplot,[0, 365],[0,0],linestyle = 2,color = 50 ;;;;;delta density plot, timemp, delta_den,ystyle = 8,xstyle = 1,ytitle='delta_Den_360km', $ yr=[-5.e-15,5.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) ;;IMF/SW ; plot, timemp, bz_intmp, ystyle = 1,xstyle = 1,ytitle='Bz',xtitle = 'DOY', $ ; yr = [-5,5], xr = [time_mp1, time_mp2],charsize = 1.5,xticks = 5 ; oplot, timemp(fp_peak),bz_intmp(fp_peak),color = 250, psym = 1, symsize = 1.5 ; plot, timemp, Vsw_intmp, ystyle = 1,xstyle = 1,ytitle='Vsw',xtitle = 'DOY', $ ; yr = [250,350], xr = [time_mp1, time_mp2],charsize = 1.5,xticks = 5 ; oplot, timemp(fp_peak),Vsw_intmp(fp_peak),color = 250, psym = 1, symsize = 1.5 ; plot, timemp, symh_intmp, ystyle = 1,xstyle = 1,ytitle='SYMH',xtitle = 'DOY', $ ; yr = [-10,10], xr = [time_mp1, time_mp2],charsize = 1.5,xticks = 5 ; oplot, timemp(fp_peak),symh_intmp(fp_peak),color = 250, psym = 1, symsize = 1.5 !p.multi = 0 device,/close set_plot,'x' ;stop maxden_n = fltarr(npeak)*!values.f_nan deltaden_n = maxden_n percent_n = maxden_n timen = maxden_n glatn = maxden_n glonn = maxden_n gltn = maxden_n mlatn = maxden_n mlonn = maxden_n mltn = maxden_n byn = maxden_n bzn = maxden_n Presn = maxden_n symhn = maxden_n maxden_s = maxden_n deltaden_s = maxden_n percent_s = maxden_n times = maxden_n glats = maxden_n glons = maxden_n glts = maxden_n mlats = maxden_n mlons = maxden_n mlts = maxden_n bys = maxden_n bzs = maxden_n Press = maxden_n symhs = maxden_n for ii = 0, npeak-1 do begin if (mlatmp(fp_peak(ii)) gt 0.0) then begin maxden_n(ii) = alog10(den_champ_360mp(fp_peak(ii))) deltaden_n(ii) = alog10(delta_den(fp_peak(ii))) percent_n(ii) = alog10(percent(fp_peak(ii))) timen(ii) = timemp(fp_peak(ii)) glatn(ii) = glatmp(fp_peak(ii)) glonn(ii) = glonmp(fp_peak(ii)) gltn(ii) = gltmp(fp_peak(ii)) mlatn(ii) = mlatmp(fp_peak(ii)) mlonn(ii) = mlonmp(fp_peak(ii)) mltn(ii) = mltmp(fp_peak(ii)) byn(ii) = by_intmp(fp_peak(ii)) bzn(ii) = bz_intmp(fp_peak(ii)) Presn(ii) = Pres_intmp(fp_peak(ii)) symhn(ii) = symh_intmp(fp_peak(ii)) endif else begin maxden_s(ii) = alog10(den_champ_360mp(fp_peak(ii))) deltaden_s(ii) = alog10(delta_den(fp_peak(ii))) percent_s(ii) = alog10(percent(fp_peak(ii))) times(ii) = timemp(fp_peak(ii)) glats(ii) = glatmp(fp_peak(ii)) glons(ii) = glonmp(fp_peak(ii)) glts(ii) = gltmp(fp_peak(ii)) mlats(ii) = mlatmp(fp_peak(ii)) mlons(ii) = mlonmp(fp_peak(ii)) mlts(ii) = mltmp(fp_peak(ii)) bys(ii) = by_intmp(fp_peak(ii)) bzs(ii) = bz_intmp(fp_peak(ii)) Press(ii) = Pres_intmp(fp_peak(ii)) symhs(ii) = symh_intmp(fp_peak(ii)) endelse endfor ;stop ;;;------------------- ;;;Write the data into a file ;;;-------------------- openw, 22,'CHAMP_360km_AACGM_peaks_30%_'+strtrim(start_end_year(jj),2)+'.dat', width = 800 printf, 22, 'time maxden delta_den percent glat glon glt mlat mlon mlt by bz Press SYMH' for i=0L, n_elements(timen)-1 do begin if (finite(timen(i))) then printf, 22, timen(i), maxden_n(i),deltaden_n(i),percent_n(i),glatn(i),glonn(i),gltn(i), $ mlatn(i),mlonn(i),mltn(i),byn(i),bzn(i), Presn(i),symhn(i) endfor for i=0L, n_elements(times)-1 do begin if (finite(times(i))) then printf, 22, times(i), maxden_s(i),deltaden_s(i),percent_s(i),glats(i),glons(i),glts(i), $ mlats(i),mlons(i),mlts(i),bys(i),bzs(i), Press(i),symhs(i) endfor close,22 ;;;set the cut for plot in polar contour fp_cut = where(maxden_n lt alog10(plot_cut)) maxden_n(fp_cut) = !Values.F_NaN fp_cut = where(maxden_s lt alog10(plot_cut)) maxden_s(fp_cut) = !Values.F_NaN ;;;;;;;;change coordinate and sort by magnitude;;--------------------- radiusn = 90.-mlatn anglen = (mltn-6.)*!pi/12. xposn = radiusn*cos(anglen) yposn = radiusn*sin(anglen) ; dcolorn = (maxden_n-alog10(abs_cut1))/(max(maxden_n,/nan)-alog10(abs_cut1))*256 dcolorn = (maxden_n-alog10(abs_cut1))/(alog10(abs_cut2)-alog10(abs_cut1))*256 fp_sortn = sort(dcolorn) xposn = xposn(fp_sortn) yposn = yposn(fp_sortn) dcolorn = dcolorn(fp_sortn) byn = byn(fp_sortn) bzn = bzn(fp_sortn) Presn = Presn(fp_sortn) radiuss = 90.+mlats angles = (mlts-6.)*!pi/12. xposs = radiuss*cos(angles) yposs = radiuss*sin(angles) ; dcolors = (maxden_s-alog10(abs_cut1))/(max(maxden_s,/nan)-alog10(abs_cut1))*256 dcolors = (maxden_s-alog10(abs_cut1))/(alog10(abs_cut2)-alog10(abs_cut1))*256 fp_sorts = sort(dcolors) xposs = xposs(fp_sorts) yposs = yposs(fp_sorts) dcolors = dcolors(fp_sorts) bys = bys(fp_sorts) bzs = bzs(fp_sorts) Press = Press(fp_sorts) endfor ;end major analysis loop. stop ;set_plot, 'ps' ;device, filename='IMF_ByBz_2005.ps',/portrait ;device,/color ;plot, by, bz,psym = 6,symsize = 0.5,xr=[-60, 60],yr=[-60, 60], $ ; xstyle =1, ystyle=1, xtitle='IMF By (nT)',ytitle='IMF Bz (nT)',title='omni 2005' ;oplot, [0,0],[-60,60] ;oplot, [-60,60],[0,0] ;;plot, by_int, bz_int,psym = 6,symsize = 0.5,xr=[-60, 60],yr=[-60, 60], $ ;; xstyle =1, ystyle=1, xtitle='IMF By (nT)',ytitle='IMF Bz (nT)',title='interpolated, 2005' ;;oplot, [0,0],[-60,60] ;;oplot, [-60,60],[0,0] ;plot, byn, bzn, psym = 6, symsize = 0.5,xr=[-60, 60],yr=[-60, 60], $ ; xstyle =1, ystyle=1, xtitle='IMF By (nT)',ytitle='IMF Bz (nT)',title='max_den_north, 2005' ;oplot, [0,0],[-60,60] ;oplot, [-60,60],[0,0] ;plot, bys, bzs, psym = 6, symsize = 0.5,xr=[-60, 60],yr=[-60, 60], $ ; xstyle =1, ystyle=1, xtitle='IMF By (nT)',ytitle='IMF Bz (nT)',title='max_den_south, 2005' ;oplot, [0,0],[-60,60] ;oplot, [-60,60],[0,0] ; ;device,/close ;set_plot,'x' ;stop set_plot,'ps' device,filename='max_CHAMP_den_360km_mp.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,'CHAMP density',charsize=1.1 xyouts,-49.0,40.0,'Northern Hemisphere',charsize=1.1 ;;;;;trajectory-------------------------------------------------------------------------- indexn = where((time ge time_mp1) and (time le time_mp2) and (mlat ge 50.),nptsn) ;NH if (nptsn gt 0) then begin mlatn_orbit = mlat(indexn) mltn_orbit = mlt(indexn) for kk = 0L, nptsn-2 do begin ;Take care of MLT time singularity. if (abs(mltn_orbit(kk)-mltn_orbit(kk+1)) gt 12) then begin ikk = where(mltn_orbit ge 12) mltn_orbit(ikk) = mltn_orbit(ikk)-24. endif endfor ineg = where(mltn_orbit lt 0.,nptneg) ;Make MLT positive if(nptneg gt 0) then mltn_orbit(ineg) = mltn_orbit(ineg) + 24. ;loadct, 33 ;plots,(90.-mlatn_orbit)*cos((mltn_orbit-6.)*!pi/12.),(90.-mlatn_orbit)*sin((mltn_orbit-6.)*!pi/12.),color = 0,psym = 7, symsize = 0.3 loadct, 0 plots,(90.-mlatn_orbit)*cos((mltn_orbit-6.)*!pi/12.),(90.-mlatn_orbit)*sin((mltn_orbit-6.)*!pi/12.),color = 210,psym = 7, symsize = 0.5 ; ;;FPI @ Resolute Bay: 74.73N, -94.89E ; glat_fpi = 74.73 ; glon_fpi = -94.89 ; height_fpi = 250. ;km ;;;;;;;;convert geographic coor to AACGM coor;;;;;;;;;;;;;;;;;;;;; ;;IDL>.com aacgmidl.pro ;;IDL>.com cnvtime.pro ;;IDL>aacgmidl ;;IDL>load_coef,"./AACGM_coef/aacgm_coeffs2010.asc" ; mlat_aacgm_fpi = fltarr(nptsn)*!Values.F_NaN ; mlon_aacgm_fpi = fltarr(nptsn)*!Values.F_NaN ; mlt_aacgm_fpi = fltarr(nptsn)*!Values.F_NaN ; ; year_fpi = 2005.;fix(date/1.e3) ; yrsec = time(indexn)*3600*24. ;sec past the start of the year ; for ia = 0, nptsn-1 do begin ; cnv_aacgm,glat_fpi, glon_fpi,height_fpi,mlat_out, mlon_out ; mlat_aacgm_fpi(ia) = mlat_out ; mlon_aacgm_fpi(ia) = mlon_out ; mlt_aacgm_fpi(ia) = calc_mlt(year_fpi, yrsec(ia), mlon_out) ; endfor ; loadct, 39 ; plots,(90.-mlat_aacgm_fpi)*cos((mlt_aacgm_fpi-6.)*!pi/12.),(90.-mlat_aacgm_fpi)*sin((mlt_aacgm_fpi-6.)*!pi/12.),color = 50,psym = -2, symsize = 0.5 endif 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. loadct, 33 index_finite = where((finite(dcolorn) eq 1),nfinite) ;index_finite = where((finite(dcolorn) eq 1) and (byn lt -10.) and (bzn lt 0.0),nfinite) xpos = xposn(index_finite) ypos = yposn(index_finite) dcolor = dcolorn(index_finite) for in = 0, nfinite-1 do begin usersym,[-1,1,1,-1,-1],[1,1,-1,-1,1],/fill,color = dcolor(in);filled square if (sqrt(xpos(in)*xpos(in)+ypos(in)*ypos(in)) le 40.) then begin plots, xpos(in), ypos(in), psym = 8, symsize = 0.8 endif endfor ;;colorbar number_levels = 33 levels_contour = dblarr(number_levels) colors_contour = dblarr(number_levels) color_table_offset =.0 maxData = alog10(abs_cut2);max(maxden_n,/nan);1.5e-15;(Aug) ;200.;max(binave_Ti); max(alog10(temps),/nan) minData = alog10(abs_cut1);1.5e-16;(Aug) ; 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, [[10^levels_contour],[10^levels_contour]],10^levels_contour,indgen(2),charsize=1.1,charthick=1, $ levels=10^levels_contour,c_color=colors_contour,/fill,position=[0.55,0.05,0.85,0.08],/noerase,ystyle=4, $ xr=[10^minData,10^maxData],xstyle=1,xticks = 2,/xlog loadct, 39 polar_dial_s,1,2 ;xyouts,.0,54.0,title1,alignment=0.5,charsize=1.2 xyouts,-49.0,44.0,'CHAMP density',charsize=1.1 xyouts,-49.0,40.0,'Southern Hemisphere',charsize=1.1 ;;;;;trajectory indexs = where((time ge time_mp1) and (time le time_mp2) and (mlat le -50.),nptss) ;SH if (nptss gt 0) then begin mlats_orbit = mlat(indexs) mlts_orbit = mlt(indexs) for kk = 0L, nptss-2 do begin ;Take care of MLT time singularity. if (abs(mlts_orbit(kk)-mlts_orbit(kk+1)) gt 12) then begin ikk = where(mlts_orbit ge 12) mlts_orbit(ikk) = mlts_orbit(ikk)-24. endif endfor ineg = where(mlts_orbit lt 0.,nptneg) ;Make MLT positive if(nptneg gt 0) then mlts_orbit(ineg) = mlts_orbit(ineg) + 24. ;;loadct, 33 ;;plots,(90.+mlats_orbit)*cos((mlts_orbit-6.)*!pi/12.),(90.+mlats_orbit)*sin((mlts_orbit-6.)*!pi/12.),color = 0,psym = 7, symsize = 0.3 loadct, 0 plots,(90.+mlats_orbit)*cos((mlts_orbit-6.)*!pi/12.),(90.+mlats_orbit)*sin((mlts_orbit-6.)*!pi/12.),color = 210,psym = 7, symsize = 0.5 ; ;;FPI @ Resolute Bay: 74.73N, -94.89E ; glat_fpi = 74.73 ; glon_fpi = -94.89 ; height_fpi = 250. ;km ;;;;;;;;convert geographic coor to AACGM coor;;;;;;;;;;;;;;;;;;;;; ;;IDL>.com aacgmidl.pro ;;IDL>.com cnvtime.pro ;;IDL>aacgmidl ;;IDL>load_coef,"./AACGM_coef/aacgm_coeffs2010.asc" ; mlat_aacgm_fpi = fltarr(nptss)*!Values.F_NaN ; mlon_aacgm_fpi = fltarr(nptss)*!Values.F_NaN ; mlt_aacgm_fpi = fltarr(nptss)*!Values.F_NaN ; ; year_fpi = 2005.;fix(date/1.e3) ; yrsec = time(indexs)*3600*24. ;sec past the start of the year ; for ia = 0, nptss-1 do begin ; cnv_aacgm,glat_fpi, glon_fpi,height_fpi,mlat_out, mlon_out ; mlat_aacgm_fpi(ia) = mlat_out ; mlon_aacgm_fpi(ia) = mlon_out ; mlt_aacgm_fpi(ia) = calc_mlt(year_fpi, yrsec(ia), mlon_out) ; endfor ;loadct, 39 ; plots,(90.-mlat_aacgm_fpi)*cos((mlt_aacgm_fpi-6.)*!pi/12.),(90.-mlat_aacgm_fpi)*sin((mlt_aacgm_fpi-6.)*!pi/12.),color = 50,psym = -2, symsize = 0.5 endif 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. loadct, 33 index_finite = where((finite(dcolors) eq 1),nfinite) ;index_finite = where((finite(dcolors) eq 1) and (bys lt -10.) and (bzs lt 0.0),nfinite) xpos = xposs(index_finite) ypos = yposs(index_finite) dcolor = dcolors(index_finite) for in = 0, nfinite-1 do begin usersym,[-1,1,1,-1,-1],[1,1,-1,-1,1],/fill,color = dcolor(in);filled square if (sqrt(xpos(in)*xpos(in)+ypos(in)*ypos(in)) le 40.) then begin plots, xpos(in), ypos(in), psym = 8, symsize = 0.8 endif endfor ;;colorbar number_levels = 33 levels_contour = dblarr(number_levels) colors_contour = dblarr(number_levels) color_table_offset =.0 maxData = alog10(abs_cut2);max(maxden_s,/nan);1.5e-15;(Aug) ;200.;max(binave_Ti); max(alog10(temps),/nan) minData = alog10(abs_cut1);1.5e-16;(Aug) ; 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, [[10^levels_contour],[10^levels_contour]],10^levels_contour,indgen(2),charsize=1.1,charthick=1, $ levels=10^levels_contour,c_color=colors_contour,/fill,position=[0.55,0.05,0.85,0.08],/noerase,ystyle=4, $ xr=[10^minData,10^maxData],xstyle=1,xticks = 2,/xlog device,/close set_plot,'x' end