vectplt.pro (Documentation for /home/davin/idl/socware/)

;$Author: nikos $
;$Date: 2014-09-03 15:05:59 -0700 (Wed, 03 Sep 2014) $
;$Header: /home/cdaweb/dev/control/RCS/vectplt.pro,v 1.16 2012/05/18 21:33:37 johnson Exp johnson $
;$Locker: johnson $
;$Revision: 15739 $
;+
; 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