;$Author: nikos $ ;$Date: 2021-02-21 18:49:50 -0800 (Sun, 21 Feb 2021) $ ;$Header: /home/cdaweb/dev/control/RCS/vectplt.pro,v 1.16 2012/05/18 21:33:37 johnson Exp johnson $ ;$Locker: johnson $ ;$Revision: 29691 $ ;+ ; NAME: VECTPLT.PRO ; ; PURPOSE: Multi-purpose vector plotting routine for CDF data file ; representation ; ; CALLING SEQUENCE: ; ; vectplt,mlats,mlons,malts,vest,vnrt,mytimes,qflgs,Stitle=station,$ ; mcors=mltin,Qmin=qmin,Qmax=qmax,nopolar=nop,Scale=scale,nobin=nob,$ ; p0lat=p0lat,p0lon=p0lon,rot=rot,binlat=bin_lat,binlon=bin_lon,$ ; limit=limit,latdel=latdel,londel=londel,Alt=alt,Ttitle=thetitle,$ ; lthik=lthik,symsiz=symsiz,symcol=symcol,_extra=extras ; ; VARIABLES: ; ; Input: ; ; mlats(*) - an N element array of geographic latitudes ; mlons(*) - an N element array of geographic longitudes ; malts(*) _ an N element array of geographic altitudes ; vest(*) - an N element array of the eastward component velocity vector ; vnrt(*) - an N element array of the northward component velocity vector ; mytimes - an IDL structure comprised of the following: year, day of year, ; month, day of month, and an N element array of seconds of day ; qflgs - an N element array of the quality flag ; ; Keyword Parameters: ; ; Stitle=station - Observing Station Name ; mcors=mltin - Coordinate Transformation Flag ; 0 - Eccentric Dipole Magnetic Local Time ; 1 - Altitude Adjusted Corrected Geomagnetic ; Coordinates (AACGC) MLT ; 2 - Geographic Coordinates ; Qmin=qmin - Minimum exceptable quality flag ; Qmax=qmax - Maximum exceptable quality flag ; nopolar=nop - Disables clockdial display; enables user defined ; projections ; Scale=scale - Set the size of vectors plotted (1000 default) ; nobin=nob - turns off binning and averaging of vectors ; binlat=bin_lat - latitude bin interval ; binlon=bin_lon - longitude bin interval ; p0lat=p0lat - map_set argument; latitude center of map ; p0lon=p0lon - map_set argument; longitude center of map ; rot=rot - map_set argument; rotation of map ; limit=limit - map_set limits ; latdel=latdel - latitude interval ; londel=londel - longitude interval ; Alt=alt - Altitude of coordinate transformation ; Ttitle=thetitle - Title for plot ; lthik=lthik - vector line thickness ; symsiz=symsiz - vector position marker ; symcol=symcol - vector position marker color ; ; ; REQUIRED PROCEDURES: ; This procedure requires the routine GENLIB.PRO ; and a shared object module, LIB_PGM.so, of C and fortran source. ; ;------------------------------------------------------------------- ; History ; ; Revision 1.1 94/10/31 09:49:24 09:49:24 baker (Kile Baker S1G) ; ; Initial version: Based on dn_cdf.pro ; ; 1.0 R. Baldwin HSTX 9/28/95 ; Initial version: Based on dn_cdf.pro 10/31/94 (Kile Baker) ; 1.1 R. Baldwin HSTX 10/23/95 ; Options included for eccentric & AACGC MLT coordinate ; transformation; geographic coordinates; polar and non-polar ; representations; vector magnitude; quality flag bounds; ; option for binning; ; 1.2 R. Baldwin HSTX 11/15/95 ; Title, margin options added; londel & latdel variables ; revised ; 1.3 R. Baldwin HSTX 12/18/95 ; Rotation of vectors is a magnetic system from ; geographic to geomagnetic ; ;------------------------------------------------------------------- function vectplt,mlats,mlons,malts,vest,vnrt,mytimes,qflgs,Stitle=station,$ mcors=mltin,Qmin=qmin,Qmax=qmax,nopolar=nop,Scale=scale,nobin=nob,$ p0lat=p0lat,p0lon=p0lon,rot=rot,binlat=bin_lat,binlon=bin_lon,$ limit=limit,latdel=latdel,londel=londel,Alt=alt,xmargin=xmargin,$ ymargin=ymargin,Ttitle=thetitle,lthik=lthik,symsiz=symsiz,symcol=symcol,$ _extra=extras ; Establish error handler catch, error_status if(error_status ne 0) then begin if((error_status eq -108) or (error_status eq -133)) then begin print, 'STATUS= No plottable data.' print, 'ERROR=Error number: ',error_status,' in vectplt.' print, 'ERROR=Error Message: ', !ERR_STRING endif else begin print, 'STATUS= Data cannot be plotted.' print, 'ERROR=Error number: ',error_status,' in vectplt.' print, 'ERROR=Error Message: ', !ERR_STRING endelse return,-1 endif icnt=n_elements(mlats) if (n_elements(mlons) ne icnt or n_elements(malts) ne icnt or $ n_elements(vest) ne icnt or n_elements(vnrt) ne icnt or $ n_elements(mytimes.times) ne icnt or n_elements(qflgs) ne icnt) then $ message, 'Arrays have incorrect size.' if(n_elements(nop) eq 0) then nop=0 if(n_elements(lthik) eq 0)then lthik = 1.7 if(n_elements(symcol) eq 0)then symcol = 100 if(n_elements(symsiz) eq 0)then symsiz = 0.5 ; Set Scale if(n_elements(scale) ne 0) then Scale=scale else scale=1000 ; Set Title yrst=strmid(strcompress(string(mytimes.year),/remove_all),2,2) monst=string(mytimes.mon) if mytimes.mon lt 10 then monst = '0'+monst dayst=string(mytimes.day) if mytimes.day lt 10 then dayst = '0'+dayst ; patch if(!d.name eq 'X') then bkgrd=255 & dtxt=0.5 if(!d.name eq 'Z') then bkgrd=0 & dtxt=1.0 ; new code if n_elements(mltin) eq 0 then mltin=0 if mltin eq 2 then nop=1 if(keyword_set(qmin)) then Qmin=qmin else qmin=0 if(keyword_set(qmax)) then Qmax=qmax else qmax=4 if(n_elements(limit) ne 0) then begin ymin=limit[0] & xmin=limit[1] & ylim=limit[2] & xlim=limit[3] endif else begin ymin=60. & xmin=-180. & ylim=90. & xlim=180. endelse ; Check for S polar plot spole=0 if(ymin eq -90. and nop eq 0) then spole=1 if(n_elements(latdel) eq 0) then latdel=10.0 if(n_elements(londel) eq 0) then londel=45.0 ; Set Binning intervals if keyword_set(nobin) then begin if(n_elements(bin_lat) ne 0) then bin_lat=bin_lat else bin_lat=1. if(n_elements(bin_lon) ne 0) then bin_lon=bin_lon else bin_lon=2. endif ; Set Altitude if(n_elements(alt) ne 0) then Alt=alt else alt=400. ;if(mltin eq 1) then begin ;for i=0,icnt-1 do begin ; malts(i)=alt ;opos=cnvcoord(mlats(i),mlons(i),malts(i)) ;mlats(i) = opos(0) ;mlons(i) = opos(1) ;endfor ;endif mglons=mlons vcount=0 if mltin eq 1 then begin print, 'AACGC ', mltin isoy=0L FOR I=0,ICNT-1 DO BEGIN malts[i]=alt opos=cnvcoord(mlats[i],mlons[i],malts[i]) mlats[i] = opos[0] mlons[i] = opos[1] mglons[i] = opos[1] temp = mytimes.times[i]+(mytimes.doy-1)*24*3600 isoy = long(temp) if( spole eq 0) then begin mlons[i] = 180.0 + 15.*mlt(mytimes.year,isoy,mlons[i]) endif else begin mlons[i] = 15.*mlt(mytimes.year,isoy,mlons[i]) endelse if mlons[i] gt 180. then mlons[i]=mlons[i]-360. ; printf,4, mytimes.year,isoy,mytimes.times(i),mytimes.doy,mlons(i) endfor utalon = (mytimes.times[0]/3600.0)*15.0 zlon = mlons[0] - utalon print, zlon,'0 UT',utalon noonstr='0 UT' endif if mltin eq 0 then begin ; print, 'ECC ', mltin sod=0L for i=0,icnt-1 do begin malts[i]=malts[i]+6371.2 sod = long(mytimes.times[i]) ;printf, 4, mytimes.year,mytimes.doy,sod,malts(i),mlats(i),mlons(i),vest(i),$ ; vnrt(i),qflgs(i) opos = eccmlt(mytimes.year,mytimes.doy,sod,malts[i],mlats[i],mlons[i]) mglons[i] = opos[0] mlats[i] = opos[1] mlons[i] = opos[2] if( spole eq 0) then begin mlons[i]= 180.0 + mlons[i]*15.0 endif else begin mlons[i]= mlons[i]*15.0 endelse iF mlons[i] gt 180. then mlons[i]=mlons[i]-360. endfor utalon = (mytimes.times(0)/3600.0)*15.0 zlon = mlons[0] - utalon ; print, zlon,'0 UT',utalon noonstr='0 UT' endif if(keyword_set(nob)) then begin b_vest=vest b_vnrt=vnrt angs = b_vest lats = mlats lons = mlons b_cont=vnrt b_cont[*]=1 b_data=b_cont b_data[*]=0 nvects = icnt b_mlat=mlats b_mlon=mglons b_malt=malts endif else begin bin_lat = 1. bin_lon = 2. nbin_lat = abs(ylim-ymin)/bin_lat nbin_lon = abs(xlim-xmin)/bin_lon lats = fltarr(nbin_lat*nbin_lon) lons = fltarr(nbin_lat*nbin_lon) ; print, ylim,ymin,xlim,xmin,nbin_lat,nbin_lon ; print, nop,mltin,qmin,qmax for j = 0,nbin_lat-1 do begin lat = float(j)*bin_lat + ymin + bin_lat/2. for i=0,nbin_lon-1 do begin lon = float(i)*bin_lon + bin_lon/2. lats[i+j*nbin_lon] = lat lons[i+j*nbin_lon] = lon - 180. endfor endfor nvects = nbin_lat*nbin_lon b_vest = fltarr(nvects) b_vnrt = fltarr(nvects) b_cont = lonarr(nvects) b_data = lonarr(nvects) b_mlat = fltarr(nvects) b_mlon = fltarr(nvects) b_malt = fltarr(nvects) b_rot = fltarr(nvects) ;patch ; print, ymin, bin_lat, bin_lon, nbin_lon for i = 0,icnt-1 do begin lat_bin = abs(fix((mlats[i]-ymin)/bin_lat)) lon_bin = fix((mlons[i]+180.)/bin_lon) indx = fix(lon_bin + lat_bin*nbin_lon) if qflgs[i] ge qmin and qflgs[i] le qmax then begin ; printf,4, indx,i,b_vest[indx],vest(i) b_vest[indx] = b_vest[indx] + vest[i] b_vnrt[indx] = b_vnrt[indx] + vnrt[i] b_cont[indx] = b_cont[indx] + 1 b_mlat[indx] = mlats[i] b_mlon[indx] = mglons[i] b_malt[indx] = malts[i] endif if qflgs[i] gt qmax then b_data[indx]=1 endfor angs = fltarr(nvects) w = where(b_cont ne 0,wc) if( wc gt 0) then begin b_vest[w] = b_vest[w]/b_cont[w] b_vnrt[w] = b_vnrt[w]/b_cont[w] endif endelse ; Compute angles w = where(b_vnrt ne 0,wc) angs[w] = atan(b_vest[w],b_vnrt[w]) ; Compute vector rotation angles for desired coordinate system ; print, 'Computing angle adjustment' for i=0, nvects-1 do begin if(b_cont[i] ne 0) then begin b_rot[i]=angadj(mltin,mytimes.year,b_mlat[i],b_mlon[i],b_malt[i]) b_rot[i]=b_rot[i]*(3.141592/180.) endif endfor w=where(b_rot ne 0, wc) if(wc gt 0) then angs[w]=angs[w]-b_rot[w] ; Compute magnitudes mags = sqrt(b_vest^2 + b_vnrt^2) if keyword_set(nop) then begin map_set,p0lat,p0lon,rot,limit=[ymin,xmin,ylim,xlim],/GRID,$ xmargin=xmargin,ymargin=ymargin,glinethick=0.5,/cyl,$ color=bkgrd,$ /noborder,latdel=latdel,londel=londel,_extra=extras if(mltin eq 2) then begin axis,xmin,ymin,xax=0,xrange=[xmin,xlim],xsty=1 axis,xmin,ymin,yax=0,yrange=[ymin,ylim],ysty=1 endif else begin axis,xmin,ymin,xax=0,xrange=[0,24],xsty=1 axis,xmin,ymin,yax=0,yrange=[ymin,ylim],ysty=1 endelse endif else begin map_set,p0lat,p0lon,rot,limit=[ymin,xmin,ylim,xlim],/GRID,$ xmargin=xmargin,ymargin=ymargin,glinethick=0.5,/stereo,$ color=bkgrd,$ /noborder,latdel=latdel,londel=londel,_extra=extras ; MLT scale for j=0,3 do begin tstr = strtrim(string(6*j),2) if(spole eq 0) then begin lon=-180+j*90 xyouts,lon,ymin-2,tstr,alignment=0.5 endif else begin lon=0-j*90 xyouts,lon,ylim+2,tstr,alignment=0.5 endelse endfor ; Invariant Latitude scale for j=0,3 do begin lat= fix(ymin) + 10*j tstr = strtrim(string(lat),2) lon=zlon xyouts,lon,lat,tstr,alignment=0.5 endfor ; UT scale for j=0,7 do begin utstr=strtrim(string(3*j),2) lon=zlon+j*45.0 if(spole eq 0) then lat=78. else lat=-78. xyouts,lon,lat,utstr,alignment=0.5 endfor ; Make tick marks on both UT and MLT scales for j=0,23 do begin lon=zlon+j*15. olon=j*15 if(spole eq 0) then begin plots,[lon,lon],[79.5,80.5] plots,[olon,olon],[ymin,(ymin+0.5)] endif else begin plots,[lon,lon],[-79.5,-80.5] plots,[olon,olon],[ylim,(ylim-0.5)] endelse endfor ; Axis if(spole eq 0) then begin utlat=replicate(80.0,360) endif else begin utlat=replicate(-80.0,360) endelse utlon=zlon+findgen(360.0) oplot,utlon,utlat,line=0 ; Axis plots,[zlon,zlon],[ymin,ylim],thick=1.7 ; Labels xyouts,0.48,0.07,'MLT',/normal xyouts,0.50,0.53,'UT',/normal endelse ; vectors determined S = findgen(16)*(!PI*2/16.) usersym,cos(S),sin(S),/FILL Re = 6362. clon = (xlim+xmin)/2. clat = (ylim+ymin)/2. side_scale = 120./scale cside = (90. - clat)*3.141592/180. bside = scale/Re*side_scale ang = 90.*3.141592/180. arg = cos(bside)*cos(cside) + sin(bside)*sin(cside)*cos(ang) if(arg lt -1)then arg=-1. if(arg gt 1)then arg=1. aside = acos(arg) tlat = 90. - aside*180./3.141592 arg = (cos(bside)-cos(aside)*cos(cside))/(sin(aside)*sin(cside)) if(arg lt -1)then arg=-1. if(arg gt 1)then arg=1. bang = acos(arg)*180./3.141592 tlon = clon + bang nx = convert_coord(clon,clat,/data,/to_normal) ndx = convert_coord(tlon,tlat,/data,/to_normal) delta = sqrt((ndx[0]-nx[0])^2+(ndx[1]-nx[1])^2) side_scale = side_scale*.08/delta cside = (90. - lats)*3.141592/180. bside = mags/Re*side_scale arg = cos(bside)*cos(cside) + sin(bside)*sin(cside)*cos(angs) arg = arg < 1. > (-1.) aside = acos(arg) tlat = 90. - aside*180./3.141592 arg = (cos(bside)-cos(aside)*cos(cside))/(sin(aside)*sin(cside)) arg = arg < 1. > (-1.) bang = acos(arg)*180./3.141592 q = where (angs lt 0,icnt) if icnt ne 0 then bang[q] = -1.*bang[q] tlon = lons + bang for i = 0,nvects-1 do begin if lons[i] gt xmin and lons[i] lt xlim and $ lats[i] gt ymin and lats[i] lt ylim and $ b_cont[i] ne 0 then begin ; ind = max(where(LVL le mags[i])) + 1 plots,lons[i],lats[i],psym=8,symsize=symsiz,color=symcol plots,[lons[i],tlon[i]],[lats[i],tlat[i]]$ ,thick=lthik endif if lons[i] gt xmin and lons[i] lt xlim and $ lats[i] gt ymin and lats[i] lt ylim and $ b_data[i] ne 0 then begin plots,lons[i],lats[i],psym=3 endif endfor plots,.84,.92,psym=8,symsize=symsiz,/normal,color=symcol plots,[.84,.92],[.92,.92],thick=lthik,/normal scstring = strtrim(string(scale),2)+'m/s' xyouts,.84,.88,scstring,/normal if(thetitle eq ' ') then begin if mltin eq 0 then begin thetitle='Velocity vs Eccentric MLT/UT and Magnetic Latitude' endif if mltin eq 1 then begin thetitle='Velocity vs AACGC MLT/UT and Magnetic Latitude' endif if mltin eq 2 then begin thetitle='Geographic Velocity ' endif endif xyouts,0.15,0.95,thetitle,/normal if(spole eq 1) then station=strmid(station,0,14) xyouts,0.05,0.85,station,/normal date = strcompress(string(fix(mytimes.mon))+'/'$' +string(fix(mytimes.day))+$ '/'+string(fix(mytimes.year)),/remove_all) xyouts,0.05, 0.9,date,/normal time_string=systime() disclaimer="Key Parameter and Survey data (labels K0,K1,K2) are preliminary data." disclaimer1="Generated by CDAWeb on: "+time_string xyouts,0.01,0.035,disclaimer,charsize=dtxt,/normal xyouts,0.01,0.01,disclaimer1,charsize=dtxt,/normal end