; ================================================================== ; Program OM_polar.pro reads MAGSAT or Oersted ASCII data file and plot the data ; "day-by-day" in the "CGM Latitude - MLT" or "Sunward - Dusk-to-Dawnward" coordinates ; (X+ magnetic north, Y+ magnetic east, Z+ downward, to the Earth's center) ; (SW+ sunward, DW+ from dusk to dawn, Z+ downward, to the Earth's center) ; The cross-polar "vector" plot comes out exactly the same in both coordinate ; systems; the time series at the bottom of the plot differ. ; Written by V. Papitashvili and F. Christiansen, DMI, April 4-5, 1999; modified ; for Oersted data on April 14, 1999 s_date=99999999 openr,42,'dato' readf,42, s_date close,42 ;s_date = 19990501 e_date = s_date ;e_date = 19990425 s_time = 000000 e_time = 235959 nrun = 0 ; should be zero for the furst run of IZMEM contours nhem = 0 ; Passes over 0 - both, 1 - north, 2 - south polar caps nfit = 0 ; 0 - no fit, 1 - linear fit to data below 60 CGM lat. nsystem = 0 ; 0 for unix, 1 for PC ndevice = 1 ; 0 for screen, 1 for PS ncoor = 1 ; 0 for CGM, 1 for SDW nwrite = 0 ; 0 for no action, 1 for writing data (only single pass) nmlt = 9 ; 9 for dependent on season, 10 for independent nangle = 12 ; 12 for oval, 11 for azimuthal nboxcar = 3 ; size of smoothing window narrow = 8 ; every Narrow point is plotted nscale = 20 ; 20 size of the vectors in cross polar plot nvect = 0 ; 0 - suppress coordinate orientation, 1 - plot nplot = 60 ; 45 or 60 degrees polar plots nmod = 0 ; 0 - mo modeled FACs, 1 - plot FACs from IZMEM/DMSP ; sat = 'MAGSAT' sat = 'oersted' !P.Font=-1 !P.Charsize=1.8 !P.color=0 !P.background=1 !P.multi = [0,0,5] oldhem='S' date = string(format='(i8)',s_date) old_date = strmid(date,0,4) + '_' + strmid(date,4,2) + '_' + strmid(date,6,2) date = old_date reads, date, format = '(i4,1x,i2,1x,i2)', iyr, imo, idy day_m = [31,28,31,30,31,30,31,31,30,31,30,31] if (iyr eq 4*floor(iyr/4)) then day_m(1) = 29 if(nsystem eq 0) then begin ; path = '/home/vp/magsat/' path='/home/fch/mag_pert/oersted/ascfiles/' endif else begin path = 'H:\mag_pert\oersted\ascfiles\' ; PC ; path = '' endelse inpf = path + date+'.asc' if(ncoor eq 0) then begin outps = date + '_cgm.ps' x_tit = 'X, nT' y_tit = 'Y, nT' endif else begin outps = date + '_sdw2.ps' x_tit = 'Sw, nT' y_tit = 'Dw, nT' endelse z_tit = 'Z, nT' if(ndevice eq 0) then begin set_plot, 'X' window,0,xsize=400,ysize=566 endif else begin set_plot, 'PS' ; device, /color device, xsize = 17., ysize = 25., xoffset = 2., yoffset = 2.5, $ /portrait, bits_per_pixel=8, filename = outps endelse ;loadct, 13 nrec = 2000 beg_file = 0 end_file = 0 openr, 1, inpf print, 'Current date: ', date tdat_row = fltarr(6) data_row = fltarr(13) time = dblarr(4,nrec) data = fltarr(13,nrec) dfit = fltarr(4,nrec) rd_newf: end_int=0 err=0 new_date = '' syr = strmid(date,0,4) smo = strmid(date,5,2) sdy = strmid(date,7,2) rd_curf: ; Read only the first row in the very first file if (beg_file eq 0) then begin ; readf, 1, format = '(i4,4i3,f10.3,3f8.3,3f8.1,3f8.3,2f7.3,2f9.3)', tdat_row,data_row readf, 1, tdat_row,data_row c_date = long(tdat_row(0)*10000)+long(tdat_row(1)*100)+long(tdat_row(2)) if (c_date lt s_date) then goto, rd_curf c_time = tdat_row(3)*10000.d0 + tdat_row(4)*100.d0 + tdat_row(5) if (c_time lt s_time) then goto, rd_curf beg_file = 1 endif data(*,0) = data_row time(1:3,0) = tdat_row(3:5) time(0,0) = tdat_row(3)*3600 + tdat_row(4)*60 + tdat_row(5) j = 0 jf = 0 rd_inpf: ON_IOERROR, endfile ;readf, 1, format = '(i4,5i3,3f8.3,3f8.1,3f8.3,2f7.3,2f9.3)', tdat_row, data_row readf, 1, tdat_row,data_row c_date = long(tdat_row(0)*10000)+long(tdat_row(1)*100)+long(tdat_row(2)) c_time = tdat_row(3)*10000.d0 + tdat_row(4)*100.d0 + tdat_row(5) ; Check end time if(c_time gt e_time) then end_int=1 ; Check if AACGM_Lat is still in the same hemisphere (above |45| degrees) abs_lat = abs(data_row(8)-data(8,j)) if (abs_lat gt 45.) then goto, continue time_cur = tdat_row(3)*3600 + tdat_row(4)*60 + tdat_row(5) if (time_cur lt time(0,j)) then time_cur = time_cur + 86400 time_dif = time_cur - time(0,j) if (abs_lat lt 45. and time_dif gt 1800) then goto, continue j = j + 1 data(*,j) = data_row time(1:3,j) = tdat_row(3:5) time(0,j) = tdat_row(3)*3600 + tdat_row(4)*60 + tdat_row(5) ; Add seconds to the current day if a part of the orbit is taken from the next day if (end_file eq 1) then time(0,j) = time(0,j) + 86400 ; Check if there is a repeated measurement at the same second and average both if ( time(0,j) eq time (0,j-1) ) then begin time(*,j-1) = (time(*,j) + time(*,j-1))/2 data(*,j-1) = (data(*,j) + data(*,j-1))/2 j = j - 1 endif ; Collect data for the data fit below 60 deg. if (nfit eq 1 and abs(data(8,j)) le 60.) then begin dfit(0,jf) = data(3,j) dfit(1,jf) = data(4,j) dfit(2,jf) = data(5,j) dfit(3,jf) = time(0,j) jf = jf + 1 endif goto, rd_inpf ;End of input_file is reached! endfile: end_file = 1 idy = idy + 1 if (idy gt day_m(imo-1)) then begin idy = 1 imo = imo + 1 smo = string(format='(i2)',imo) if (imo lt 10) then begin smo = '0'+string(format='(i1)',imo) endif else begin smo = string(format='(i2)',imo) endelse if (imo gt 12) then begin imo = imo - 12 smo = '0'+string(format='(i1)',imo) iyr = iyr + 1 endif endif if (idy lt 10) then begin sdy = '0'+string(format='(i1)',idy) endif else begin sdy = string(format='(i2)',idy) endelse syr = string(format='(i4)',iyr) date = syr + '_' + smo + '_' + sdy new_date = ' - ' + date inpf = path + date+'.asc' close, 1 openr, 1, inpf, ERROR = err if (err NE 0) then goto, continue print, 'Start reading new date: ', date goto, rd_inpf continue: ; if ( time(0,0) eq time(0,j) ) then goto, end_prog if (nhem eq 1 and data(8,0) lt 0.) then goto, skip_plot if (nhem eq 2 and data(8,0) gt 0.) then goto, skip_plot if (j lt nboxcar) then goto, skip_plot print, 'Start',time(0,0),' and end seconds:',time(0,j),' (24 hr = 86400 s) # of points read:',j+1 ; Fitting data below 60 deg ******************************** if (nfit eq 1) then begin xbeg = dfit(3,0) for kf = 0,2 do begin ab = linfit(dfit(3,0:jf-1)-xbeg,dfit(kf,0:jf-1)) ; ab = poly_fit(dfit(3,0:jf-1)-xbeg,dfit(kf,0:jf-1),2) print, ab xcur = time(0,*) - xbeg data(kf+3,*) = data(kf+3,*) - (ab(0) + ab(1)*xcur) ; data(kf+3,*)=data(kf+3,*)-(ab(0)+ab(1)*xcur+ab(2)*xcur^2) endfor endif ; Plotting the cross-polar pass lon = fltarr(361) lat = fltarr(3,361) xr = fltarr(3,361) yr = fltarr(3,361) lon = findgen(361)*!pi/180 for i = 0,2 do begin lat(i,*) = (i + 1)*15. xr(i,*) = -lat(i,*)*sin(lon) yr(i,*) = lat(i,*)*cos(lon) endfor r = fltarr(36) ln = fltarr(36) for i = 0,35 do begin ln(i) = (float(fix(i/3))/4.)*2.*!pi if (i + 2)/3 eq (float(i) + 2.)/3. then r(i) = 15. endfor xl = 3.0*r*cos(ln) yl = 3.0*r*sin(ln) ;insert empty plot for missing hemisphere if(data(8,0) gt 0) then begin if(oldhem eq 'N') then begin plot, xr(2,*), yr(2,*), xstyle = 4, ystyle = 4, pos = [0.18, 0.58, 0.82, 1.12],$ xrange = [max(xr),min(xr)], yrange = [max(yr),min(yr)] oplot, xr(1,*), yr(1,*) oplot, xr(0,*), yr(0,*) oplot, xl, yl plot,xr(2,*),yr(2,*),xstyle=4,ystyle=4,/nodata plot,xr(2,*),yr(2,*),xstyle=4,ystyle=4,/nodata plot,xr(2,*),yr(2,*),xstyle=4,ystyle=4,/nodata plot,xr(2,*),yr(2,*),xstyle=4,ystyle=4,/nodata endif oldhem='N' endif else begin if(oldhem eq 'S') then begin plot, xr(2,*), yr(2,*), xstyle = 4, ystyle = 4, pos = [0.18, 0.58, 0.82, 1.12],$ xrange = [max(xr),min(xr)], yrange = [max(yr),min(yr)] oplot, xr(1,*), yr(1,*) oplot, xr(0,*), yr(0,*) oplot, xl, yl plot,xr(2,*),yr(2,*),xstyle=4,ystyle=4,/nodata plot,xr(2,*),yr(2,*),xstyle=4,ystyle=4,/nodata plot,xr(2,*),yr(2,*),xstyle=4,ystyle=4,/nodata plot,xr(2,*),yr(2,*),xstyle=4,ystyle=4,/nodata endif oldhem='S' endelse plot, xr(2,*), yr(2,*), xstyle = 4, ystyle = 4, pos = [0.18, 0.58, 0.82, 1.12],$ xrange = [max(xr),min(xr)], yrange = [max(yr),min(yr)] oplot, xr(1,*), yr(1,*) oplot, xr(0,*), yr(0,*) oplot, xl, yl if( time(0,j-1) lt 86400) then new_date='' xyouts, 0.0, -52.0, '!5 ' + sat +'_11a_00 ' + old_date + new_date, $ ;xyouts, 0.0, -50.0, '!5 ' + sat +' ' + old_date, $ alignment = 0.5, charsize = 1.5, /noclip old_date = date ; Different labels for North and South hemispheres splot = string(format='(i2)',nplot) if (data(8,0) gt 0.) then begin hem = '!5North' deg = '!5 '+splot+'!17_' zvec = '!20S!N' delta = 7. tm = 3. oplot, [55.,55.],[49.,49.],psym = 4 xyouts, 50., 49.5, 'ascending node', alignment = 0., charsize = 1.0, /noclip endif else begin hem = '!5South' deg = '!5-'+splot+'!17_' zvec = '!9n!N' delta = -7. tm = 15. oplot, [-55.,-55.],[49.,49.],psym = 4 xyouts, -50., 49.5, 'ascending node', alignment = 1., charsize = 1.0, /noclip endelse xyouts,-30.0, -37.0, hem, /noclip ,charsize=1.25 xyouts,-32.0, -32.0, deg, /noclip ,charsize=1.25 ;xyouts,-43.0, -15.0, hem, alignment = 0.0, /noclip ;xyouts,-45.0, -10.0, deg, alignment = 0.0, /noclip xyouts,-46.0, 7.5,'!5MLT', alignment = 0.0, /noclip ,charsize=1.25 xyouts, 0.0, 49.5, '!500', alignment = 0.5, /noclip ,charsize=1.25 xyouts,-46.0, 1.5, '!506', alignment = 0.0, /noclip ,charsize=1.25 xyouts, 0.0, -46.0, '!512', alignment = 0.5, /noclip ,charsize=1.25 xyouts, 46.0, 1.5, '!518', alignment = 1.0, /noclip ,charsize=1.25 if (nvect eq 1) then begin ; Plot X, Y, Z vectors orientation xv = fltarr(2) yv = fltarr(2) if(ncoor eq 0) then begin svec = sin(tm*15.*0.017453293) cvec = cos(tm*15.*0.017453293) xv(0) = -45. * svec yv(0) = 45. * cvec xv(1) = (-45.+delta) * svec yv(1) = ( 45.-delta) * cvec xyouts, xv(0)+1.5, yv(0)+1.5, zvec, /noclip xyouts, xv(0)-2.5, yv(0)+2.5, '!5Z', /noclip oplot, xv, yv, thick = 4 xyouts, xv(1)-1., yv(1)-1., '!5X', /noclip if (tm lt 15.) then begin xv1 = xv(0)-abs(xv(0)-xv(1)) oplot, [xv(0),xv1],yv, thick = 4 xyouts, xv1-2., yv(1)-2, '!5Y', /noclip endif else begin yv1 = yv(0)+abs(yv(0)-yv(1)) oplot, xv, [yv(0),yv1], thick = 4 xyouts, xv(1)+3., yv1+3., '!5Y', /noclip endelse endif else begin svec = sin(15.*15.*0.017453293) cvec = cos(15.*15.*0.017453293) xv(0) = -45. * svec yv(0) = 45. * cvec xv(1) = xv(0) yv(1) = yv(0) - 10. oplot, xv, yv, thick = 4 xyouts, xv(1)-1., yv(1)-1., '!5Sw', /noclip xv(1) = xv(0) - 10. yv(1) = yv(0) oplot, xv, yv, thick = 4 xyouts, xv(1)-1., yv(1)-1., '!5Dw', /noclip xyouts, xv(0)+1.5, yv(0)+1.5, zvec, /noclip xyouts, xv(0)+4., yv(0), '!5Z', /noclip endelse endif ; Convert coordinates to plot the satellite track xrm = fltarr(j) yrm = fltarr(j) sinmlt = fltarr(j) cosmlt = fltarr(j) ; Plotting the track by vectors; (diamond=4 is entry, ascending node) if (nplot eq 45) then begin npl = 1. endif else begin npl = 1.5 endelse sinmlt = sin(0.017453293*data(nmlt,*)*15.) cosmlt = cos(0.017453293*data(nmlt,*)*15.) xrm = -(90-abs(data(8,*)))*sinmlt * npl yrm = (90-abs(data(8,*)))*cosmlt * npl if (nplot eq 45) then begin oplot, [xrm(0)],[yrm(0)],psym=4 endif else begin for i = 0,j-2 do begin if(abs(data(8,i)) le float(nplot) and abs(data(8,i+1)) ge float(nplot)) then begin oplot, [xrm(i+1)],[yrm(i+1)],psym=4 goto, go_out endif endfor endelse go_out: if (nmod eq 1) then begin !P.color=80 nrun = nrun + 1 OM_cont, nhem, nrun !P.color=1 endif ; oplot, [xrm(j-1)],[yrm(j-1)],psym=2 ;for i = 1,j-2 do begin ; oplot, [xrm(i)],[yrm(i)],psym=3 ;endfor ; oplot, [-20.,-27.],[42.5,42.5],linestyle = 0 ;xyouts, -30., 43., 'classic MLT', alignment = 0., charsize = 0.8, /noclip ; for i = 0,7 do begin ; oplot, [-20.-i],[45.5+0.*i], psym=3 ; endfor ;xyouts, -30., 46., 'approx MLT', alignment = 0., charsize = 0.8, /noclip ;xyouts, -25., 45., 'approx. MLT and meridian angle', alignment = 0., charsize = 0.8, /noclip ; Now prepare to plot the data ; Defining max-min for both axes xtime = fltarr(j) datxo = fltarr(j) datyo = fltarr(j) datzo = fltarr(j) zero = fltarr(j) cosang = fltarr(j) sinang = fltarr(j) dum = fltarr(j) ; Convert data from NE to XY-geomagnetic coordinates using the Azimuth_angle (11) and ; then the Oval_angle (12). Note that data can be smoothed - e.g., by the boxcar ; window of 15 sec ; data(nangle,0:j-1) = 0. ; Used for plotting in geographic coords sinang = sin(0.017453293*data(nangle,0:j-1)) cosang = cos(0.017453293*data(nangle,0:j-1)) datxo = data(3,0:j-1)*cosang+data(4,0:j-1)*sinang dum = smooth(datxo,nboxcar) datxo = dum datyo=-data(3,0:j-1)*sinang+data(4,0:j-1)*cosang dum = smooth(datyo,nboxcar) datyo = dum datz = smooth(data(5,0:j-1),nboxcar) ; Rotate the vectors to the XY-coordinates of the plot, then scale them and add them ; to the trajectory... plot every 15th point vecxo = fltarr(j) vecyo = fltarr(j) if (data(8,0) gt 0.) then begin vecxo = -datyo*cosmlt + datxo*sinmlt vecyo = -datyo*sinmlt - datxo*cosmlt endif else begin vecxo = -datyo*cosmlt - datxo*sinmlt vecyo = -datyo*sinmlt + datxo*cosmlt endelse if(ncoor eq 1) then begin datxo=-vecyo datyo=-vecxo datz=-datz endif if(c_time gt 0 and nwrite eq 1) then begin ih=floor(c_time/10000) im=round((c_time-10000*ih)/100) if(ih lt 10) then begin strh='0'+string(format='(i1)',ih) endif else begin strh=string(format='(i2)',ih) endelse if(im lt 10) then begin strm='0'+string(format='(i1)',im) endif else begin strm=string(format='(i2)',im) endelse openw, 10,'data/'+date+':'+strh+':'+strm+'.dat' for i=0,j-1 do begin printf, 10, xrm(i),yrm(i),vecxo(i),vecyo(i) endfor close, 10 endif vecxo = vecxo/nscale+xrm vecyo = vecyo/nscale+yrm for i=1,j-2,narrow do begin if(abs(data(8,i)) ge nplot) then $ arrow, xrm(i), yrm(i), vecxo(i), vecyo(i), /data, hsize = -0. endfor ; Reverse sequence for the Southern hemisphere (plot from right to left) if (tm eq 3.) then begin xtime = time(0,0:j-1) xmin = xtime(0) xmax = xtime(j-1) endif else begin xtime = reverse(time(0,0:j-1)) xmax = xtime(0) xmin = xtime(j-1) dum = reverse(datxo) datxo = dum dum = reverse(datyo) datyo = dum dum = reverse(datz) datz = dum endelse zero(*) = 0. all = [[datxo],[datyo],[datz]] allmin = min(all) allmax = max(all) !P.charsize = 1.8 ; Dummy plot - to make a room for the dial plotted above plot, xtime, datz, /nodata, xstyle=4, ystyle =4 blank = strarr(10) blank(*) = ' ' ; Plot Magnetic North component (X) and get ticks ; (suppress "Seconds of Day" as numbers (e.g.12x10^4) and ; plot them later as characters (e.g., 12000) plot, xtime, datxo, ytitle = x_tit, psym=3, $ /xstyle, xrange=[xmin,xmax], yrange=[allmin,allmax], xtick_get = sec_day, $ xtickname = blank oplot, xtime, zero(*), linestyle = 2, thick = 1.0 ; Define extra notations for the horizontal axis nsec = size(sec_day, /n_elements) sec_d = strarr(nsec) g_lat = strarr(nsec) g_lon = strarr(nsec) a_lat = strarr(nsec) hh_mm = strarr(nsec) h_mlt = strarr(nsec) for i = 0,nsec-1 do begin hh_mm(i) = string(floor(sec_day(i)/3600),format='(i2)') + ':' + $ string(floor(sec_day(i)/60) mod 60,format='(i2)') + ':' + $ string(sec_day(i) mod 60,format='(i2)') sec_d(i) = string(sec_day(i),format='(i5)') ns = where(floor(xtime/10)*10 eq sec_day(i),ncount) if (ncount gt 0) then begin g_lat(i) = string(data(0,ns(0)),format='(f5.1)') g_lon(i) = string(data(1,ns(0)),format='(f5.1)') a_lat(i) = string(data(8,ns(0)),format='(f5.1)') h_mlt(i) = string(data(9,ns(0)),format='(f5.2)') endif else begin g_lat(i) = ' ' g_lon(i) = ' ' a_lat(i) = ' ' h_mlt(i) = ' ' endelse endfor ; Plot UT on the top and "Seconds of Day" at the bottom of the first box axis, xaxis = 1, xtickv = sec_day, xticks = nsec, xtickname = hh_mm,charsize=2. xyouts, 0.07, 0.587, 'UT', charsize = 0.9, align = 1., /normal, /noclip axis, xaxis = 0, xtickv = sec_day, xticks = nsec, xtickname = sec_d,charsize=2. xyouts, 0.07, 0.425, 'Seconds', charsize = 0.9, align = 1., /normal, /noclip ; Plot Magnetic East component (Y) and MLT (top) and AACGM latitude (bottom) plot, xtime, datyo, ytitle = y_tit, psym=3, /xstyle, $ xrange=[xmin,xmax], yrange=[allmin,allmax], xtickv=sec_day, xtickname=blank oplot, xtime, zero(*), linestyle = 2, thick = 1.0 axis, xaxis = 1, xtickv = sec_day, xticks = nsec, xtickname = h_mlt,charsize=2. xyouts, 0.07, 0.387, 'MLT', charsize = 0.9, align = 1., /normal, /noclip axis, xaxis = 0, xtickv = sec_day, xticks = nsec, xtickname = a_lat,charsize=2. xyouts, 0.07, 0.225, 'AACGM Lat', charsize = 0.9, align = 1., /normal, /noclip ; Plot Magnetic Vertical (Downward, to the Earth's center) component (Z) ; and geographic lomgitude (top) and latitude (bottom) plot, xtime, datz, ytitle = z_tit, psym=3, /xstyle, $ xrange=[xmin,xmax], yrange=[allmin,allmax], xtickv=sec_day, xtickname=blank oplot, xtime, zero(*), linestyle = 2, thick = 1.0 axis, xaxis = 1, xtickv = sec_day, xticks = nsec, xtickname = g_lon,charsize=2. xyouts, 0.07, 0.187, 'Geo Lon', charsize = 0.9, align = 1., /normal, /noclip axis, xaxis = 0, xtickv = sec_day, xticks = nsec, xtickname = g_lat,charsize=2. xyouts, 0.07, 0.025, 'Geo Lat', charsize =0.9, align = 1., /normal, /noclip skip_plot: ; Ending to print a single cross-polar orbit... go to read the next pass if (c_date gt e_date) then goto, end_prog if (end_int eq 1) then goto, end_prog ; Check where to go to continue reading the input file(s) if(ndevice eq 0) then begin read, dummy if(dummy eq 0) then goto, end_prog endif if(err ne 0) then goto,end_prog if (end_file eq 0) then begin goto, rd_curf endif ; else begin ; end_file = 0 ; goto, rd_newf ;endelse end_prog: close, 1 if(ndevice eq 1) then device, /close_file print, 'Ending...' end