;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 = 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('./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' ;;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) traj_time_75 = fltarr(ntraj+1)*!values.f_nan ;;for the first orbit traj_time = timemp(0:fp_traj(0)) traj_mlat = mlatmp(0:fp_traj(0)) fp_75 = where(abs(traj_mlat) ge 75., nmlat_75) if (nmlat_75 ge 2) then begin traj_time_75(0) = traj_time(fp_75(nmlat_75-1))-traj_time(fp_75(0)) ; time over 75 mlat, unit of day endif for itra = 0L, ntraj-2 do begin traj_time = timemp(fp_traj(itra):fp_traj(itra+1)) traj_mlat = mlatmp(fp_traj(itra):fp_traj(itra+1)) fp_75 = where(abs(traj_mlat) ge 75., nmlat_75) if (nmlat_75 ge 2) then begin traj_time_75(itra+1) = traj_time(fp_75(nmlat_75-1))-traj_time(fp_75(0)) ; time over 75 mlat endif endfor ;;for the last orbit traj_time = timemp(fp_traj(ntraj-1):nptsmp-1) traj_mlat = mlatmp(fp_traj(ntraj-1):nptsmp-1) fp_75 = where(abs(traj_mlat) ge 75., nmlat_75) if (nmlat_75 ge 2) then begin traj_time_75(ntraj) = traj_time(fp_75(nmlat_75-1))-traj_time(fp_75(0)) ; time over 75 mlat endif fp_nan = where(traj_time_75 gt 0.1) traj_time_75(fp_nan) = !values.f_nan percent_pass_75 = traj_time_75*24.*3600.*Vsc/1.e3/(2.*(360.+6371.)*(90.-75.)*!pi/180.)*100. set_plot, 'ps' device, filename='CHAMP_360km_AACGM_orbit_75MLAT'+strtrim(start_end_year(jj),2)+'.ps',/portrait device,/color plot, percent_pass_75, psym = 4,ytitle='percent of CHAMP orbit>75MLAT',xstyle=1, yr = [0, 150] xyouts, 2000., 130., 'mean = '+strmid(mean(percent_pass_75,/nan),2) device,/close set_plot,'x' stop ;;;;;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 endfor ;end major analysis loop. end