c************
c   Description: Geo primitive / plotter
c
c   Author:    JGS   8/1/90
c
c   $Revision: 1.1 $
c   $Modtime:   26 Aug 1997 12:56:18  $
c
c  Revisions:	GE Doty, 9/98  (following DN Simkins, RR Weyker, COGSX)
c		  - Include COGSX modifications
c		  - Allow rtdelta to go up to 300
c		  - Keep realtime mode working after new Zulu day
c		  - G_TIME result to pass out time-of-day
c		  - Event/cone line-width controls
c		  - /SVH for unique SV colors with matching horizons
c		  - /SVIEW, /MVIEW for event display only if site, mark in view
c		  - site- and mark-AoA spec options
c		  - workaround bug in cone plots (look for 'cbug' comment)
c		  - add 'THICK' keyword for Type 5000 non-backing store files
c		  - /XOFFS=n offsets base longitude from -180/+180 so horizon
c			lines not cut off on map edge
c
c	11-03 JLR/EHG DR 500002 For Intel Compiler, changed .not. to not and type to write.
c
c************
	program XGEOPLOT
	include 'xgeoplot.inc'

	call M$CHECKIN
	call XGEO_INIT ()

	call M$SYNC ()
	call XGEO_PLOT ()

	do while (Gx.status.ge.0 .and. .not. (Mc.break.or.Mx.breakflag))
	  Gx.status = 0
	  call XGEO_EVENT ()
	  call XGEO_MESSAGE ()
	  call XGEO_PIPES ()
	  if (Gx.status .eq. 0) call M$PAUSE(Mc.pause)
	enddo

	call XGEO_END ()
	call M$RETURN

	end

	subroutine XGEO_INIT ()
	include 'xgeoplot.inc'
	pointer (pev,ev)
c
	colors = 'Grid,Horizon,TDOA,FDOA,Trace5000,'//
     &'Coastlines,Political,InterPol,Rivers,RailRoads,Roads,UserMap'
c
c init parameters
	Gx.status = 0
	Gx.pdis = 6				! persp distance (earth radii)
	Gx.plat = 0				! persp latitude
	Gx.plon = 0				! persp longitude
	Gx.view = max(1,M$GET_SWITCH('VIEW'))   ! 1=latlon 2=persp 3=svA 4=svB
	Gx.spec = M$GET_SSWITCH_DEF('SPECS',1)  ! Display mouse specs text
	Gx.mark = M$GET_SSWITCH_DEF('MARK',1)	! Display marker
	Gx.axes = M$GET_SSWITCH_DEF('AXES',0)	! Display axes on plots
	Gx.prop = .true.			! Aspect ratio enforce
	Gx.labs = .true.			! Labels on layers
	Gx.pan  = .true.			! Enable panning on zooms
	Gx.bar  = .false.			! Disable menu bar
	Gx.hold = .false.			! Disable hold mode
	Gx.sorg = .true.			! set origin with LMB click mode
	Gx.vcon = .false.			! no vehicle cones
	Gx.zfill= .true.			! have zoom fill window
	Gx.evnt = M$GET_SSWITCH('EVENT')	! Enable events
	Gx.cep  = M$GET_SSWITCH('CEP')		! Enable cep
	Gx.split = M$GET_SSWITCH('SPLIT')	! Split screen zoom mode
	Gx.units = M$GET_SWITCH_DEF('UNITS',1)	! Metric/English/? units
	Gx.tlmax = M$GET_SWITCH_DEF('TLMAX',24)	! Text label max
	Gx.tlreg = M$GET_SWITCH_DEF('TLREG',10)	! Text label region
	i = M$TPICK(5,t1,t2)
	Gx.time = t1 + t2			! Global time of interest
	Gx.tmin = 0				! Global tmin
	Gx.tmax = 86400				! Global tmax
	Gx.zk = 1		! last zoom level
	Gx.zs = 0		! last split mode
	Gx.speed = 0.040	! min display time in seconds (25 Hz)
	Gx.glmode = 1		! geo location mode
	Gx.glevn1 = 1		! locate event 1
	Gx.glevn2 = 1		! locate event 2
	Gx.glstd1 = 1e-6	! std of measurement 1
	Gx.glstd2 = 1e-2	! std of measurement 2
	Gx.cepstep = 1000	!
	Gx.iballa = 1		! ball A
	Gx.iballb = 2		! ball B
	Gx.coneveh = -1		! all vehicles
	Gx.coneang = 1		! vehicle cone angle
	Gx.conemod = 1		! cone plotting mode
	Gx.conediv = 60		! divisions around cone ellipse
	Gx.conelev = 1		! number of cone slices
	Gx.rtdelta = M$GET_SWITCH('RT')
	Gx.rtlast  = 0
	Gx.linew   = 1				! Event/cone line width
	Gx.svh     = M$GET_SSWITCH('SVH')	! SV, horizon colors uniq, match
	Gx.sview   = M$GET_SSWITCH('SVIEW')	!Show event only if Site in view
	Gx.mview   = M$GET_SSWITCH('MVIEW')	!Show event only if Mark in view
	Gx.aoaspec = M$GET_SSWITCH('AOASPEC')	! Add site-, mark-AoA to specs
c
c init marker T5 layer
	call T5INIT (Gx.mrk, P_TEXT, PC_HORIZON, L_circlesymbol, 3)
	Gx.mrk.name_len = M$GET_USWITCH('MARKSTR',Gx.mrk.name)
c
c init event T5 layer
	do i = 1,5
	  call T5INIT (Gx.evp(i), P_TEXT, -1, L_Xsymbol, 2)
	  Gx.evp(i).name_len = 2  ! filled in by DRAW_EVENT
	enddo
c
c open X-window display and root window
	Mx.origin = 1
	call MX$GRAPHICSINIT()
	call MXF$SETUP (Mx,Mc.width,Mc.height) 
	call MX$OPEN ('X-Geoplot')
c
c rubberband for images
	call MX$DRAW_MODE (901)
c
c backing store setup
	Gx.bs = M$GET_SSWITCH_DEF('BS',1)
	do k = 1,L_stack
	  Gx.bsmap(k).ip = 0
	enddo
c
c get control mode
	Gx.cntrls = M$GET_SWITCH ('XCNT')
	if (Gx.cntrls.gt.2 .or. Gx.cntrls.lt.0) Gx.cntrls = 1
c
c define the colormap
	Gx.ncolors = M$GET_SWITCH_DEF('XCMAX',L_colors)	! number of colors
	call MX$COLORMAP (Mc.colormap(4), Gx.ncolors)
	if (M$GET_USWITCH('XCF',str).le.0) str = ' '
	call XGEO_COLORS (0,str)
c
c init the view dimensions
	Gx.xoffs = M$GET_SWITCH_DEF ('XOFFS', 0)
	call XGEO_RESET (-1)
	Gx.rlon1 = Gx.x1			! event region of interest
	Gx.rlon2 = Gx.x2
	Gx.rlat1 = Gx.y1
	Gx.rlat2 = Gx.y2
c
c init the facilities structures
	call XGEO_FACILITIES (0)
	Gx.xf.msgid   = M$GET_SWITCH ('MSGID')
	Gx.xf.msgmask = max (0, M$GET_SWITCH ('MASK'))
c
c init transformation for type 5000 files
	call XGEO_GET_DATA (0,' ',-1)
c
c init sincos lookup table for fast persp transformation
	call SINCOS_INIT ()
c
c get the site
	str = M$APICK(2)
	call XGEO_LOAD (P_SITE,str)
c
c open the state vector files
	str = M$APICK(3)
	if (M$PPICK(4)) str = str(1:M$LENGTH(str))//'|'//M$APICK(4)
	call XGEO_LOAD (P_SV,str)
c
c init the grid (1st on layer link list)
	call M$MALLOC (L_layersz,i)
	call M$MALLOC (L_grid*8,j)
	call XGEO_GET_GRID (%val(i),%val(j))
	Gx.nlayer = 0
	Gx.lloc(Gx.nlayer) = i
c
c init the map files (FILE1|FILE2|FILE3...)
	i = max(0, M$UAPICK(1,str,0) )
	call XGEO_LOAD (P_LAYER,str(1:i))
c
c init the base event (1st on event link list)
	Gx.nevent = 0
	if (Gx.nball.gt.1) then
	  Gx.nevent = 1
	  call M$MALLOC (L_eventsz,Gx.eloc(Gx.nevent)) 
	  call XGEO_EVENT_INIT (Gx.nevent)
	endif
c
c allocate dynamic memory for grids
	Gx.lmg = 0
	Gx.nmg = M$GET_SWITCH('NMG')
	if (Gx.nmg .lt. 4) Gx.nmg = L_mg
	call M$MALLOC (Gx.nmg*Gx.nmg*4*5,Gx.lmg)
c
c open the altitude database file
	if (M$GET_USWITCH('ADBF',str) .gt. 0) then
	  call XGEO_LOAD (P_ALTDB,str)
	else
	  Gx.lhadb = 0
	endif
c
c loop through the positionally independent parameters
	i = 6
	do while (M$PPICK(i))
	  str = M$UPICK(i)
	  j = M$SEARCH(str,'=')
	  if (j.eq.0) then
	    call M$WARNING ('Argument is not PosInd: '//str)	
	  elseif (M$A2D(str(j+1:),dval).le.0) then
	    call M$WARNING ('Problem parsing key value: '//str)	
	  else
	    call XGEO_SETD (str(1:j-1),0,dval)
	  endif
	  i = i+1
	enddo
c
c create results for communicating with other XGEO macros
	call M$ARSLT ('G_ACTION','DONE')
	call M$DRSLT ('G_TIME', Gx.time)
c
	return
	end



	subroutine XGEO_SETD (key,info,dval)
	include 'xgeoplot.inc'
	pointer (ply,ly)
	pointer (pev,ev)
	character key*(*)
	integer*4 info

	if     (key .eq. 'VIEW') then
	  Gx.view = dval

	elseif (key .eq. 'PAN') then
	  Gx.pan = dval .gt. 0
	elseif (key .eq. 'BS') then
	  Gx.bs = dval .gt. 0
	elseif (key .eq. 'SPLIT') then
	  Gx.split = dval .gt. 0
	elseif (key .eq. 'LABELS') then
	  Gx.labs = dval .gt. 0
	elseif (key .eq. 'AXES') then
	  Gx.axes = dval .gt. 0
	elseif (key .eq. 'ASPECT') then
	  Gx.prop = dval .gt. 0
	elseif (key .eq. 'SPECS') then
	  Gx.spec = dval .gt. 0
	elseif (key .eq. 'BAR') then
	  Gx.bar = dval .gt. 0
	elseif (key .eq. 'EVENTS') then
	  Gx.evnt = dval .gt. 0
	elseif (key .eq. 'CEP') then
	  Gx.cep = dval .gt. 0
	elseif (key .eq. 'SORG') then
	  Gx.sorg = dval .gt. 0
	elseif (key .eq. 'MARK') then
	  Gx.mark = dval .gt. 0
	elseif (key .eq. 'HOLD') then
	  Gx.hold = dval .gt. 0
	elseif (key .eq. 'CONE') then
	  Gx.vcon = dval .gt. 0
	elseif (key .eq. 'ZFILL') then
	  Gx.zfill = dval .gt. 0
	elseif (key .eq. 'UNITS') then
	  Gx.units = dval
	elseif (key .eq. 'SVH') then
	  Gx.svh = dval .gt. 0
	elseif (key .eq. 'SVIEW') then
	  Gx.sview = dval .gt. 0
	elseif (key .eq. 'MVIEW') then
	  Gx.mview = dval .gt. 0
	elseif (key .eq. 'AOASPEC') then
	  Gx.aoaspec = dval .gt. 0


	elseif (key .eq. 'LAYER' .and. info .ge. 0) then
	  lval = dval
	  ply = Gx.lloc(info)
	  if (lval.eq.0) ly.enable = JIBCLR(ly.enable,PM_GLOBAL)
	  if (lval.eq.1) ly.enable = JIBSET(ly.enable,PM_GLOBAL)
	  if (lval.gt.1) ly.enable = lval
	elseif (key .eq. 'EVENT' .and. info .gt. 0) then
	  lval = dval
	  pev = Gx.eloc(info)
	  if (lval.eq.0) ev.enable = JIBCLR(ev.enable,PM_GLOBAL)
	  if (lval.eq.1) ev.enable = JIBSET(ev.enable,PM_GLOBAL)
	  if (lval.gt.1) ev.enable = lval
	elseif (key .eq. 'XOFFS') then
	  Gx.xoffs = dval
	elseif (key .eq. 'X1') then
	  Gx.x1 = dval
	elseif (key .eq. 'X2') then
	  Gx.x2 = dval
	elseif (key .eq. 'Y1') then
	  Gx.y1 = dval
	elseif (key .eq. 'Y2') then
	  Gx.y2 = dval
	elseif (key .eq. 'EX1') then
	  Gx.rlon1 = dval
	elseif (key .eq. 'EX2') then
	  Gx.rlon2 = dval
	elseif (key .eq. 'EY1') then
	  Gx.rlat1 = dval
	elseif (key .eq. 'EY2') then
	  Gx.rlat2 = dval

	elseif (key .eq. 'RTDELTA') then
	  Gx.rtdelta = dval
	elseif (key .eq. 'CONEMOD') then
	  Gx.conemod = dval
	elseif (key .eq. 'CONELEV') then
	  Gx.conelev = dval
	  if (info.ne.0) Gx.conemod = info
	elseif (key .eq. 'CONEDIV') then
	  Gx.conediv = dval
	elseif (key .eq. 'CONEANG') then
	  Gx.coneang = dval
	elseif (key .eq. 'CONEVEH') then
	  Gx.coneveh = dval
	elseif (key .eq. 'TLREG') then
	  Gx.tlreg = dval
	elseif (key .eq. 'TLMAX') then
	  Gx.tlmax = dval
	elseif (key .eq. 'MARKSYMB') then
	  Gx.mrk.symbol = dval
	elseif (key .eq. 'MARKSIZE') then
	  Gx.mrk.size = dval
	elseif (key .eq. 'MARKCOLO') then
	  Gx.mrk.color = dval
	elseif (key .eq. 'NMG') then
	  Gx.nmg = max(3.d0,dval)
	  call M$REMALLOC (Gx.nmg*Gx.nmg*4*5,Gx.lmg)
	else
	  call M$WARNING ('Unrecognized key: '//key)
	endif
	return
	end


	subroutine XGEO_LOAD (type,string)
	include 'xgeoplot.inc'
	pointer (ply,ly)
	record /HEADER/ hcb
	integer*4 type
	character string*(*)

	if (type.eq.P_SITE) then			! site
	  call M$INIT (hcb,string,'5000',' ',0)
	  call M$OPEN (hcb,HCBF_INPUT+HCBF_OPTIONAL)
	  if (hcb.open) call M$GRAB5N (hcb,Gx.site,1.d0,1,-1,'POS')
	  call M$CLOSE (hcb)

	elseif (type.eq.P_SV) then			! state vectors
	  str = string
	  Gx.nball = 0
	  do while (M$LENGTH(str).gt.0)
	    j = M$SEARCH(str,'|')
	    if (j.eq.0) j = M$LENGTH(str)+1
	    hcb.file_name = str(1:j-1)
	    i = M$APICK_POST (hcb.file_name,j-1)
	    call M$INIT (hcb,hcb.file_name,'5000',' ',0)
	    call M$OPEN (hcb,HCBF_INPUT+HCBF_NOABORT)
	    if (.not. hcb.open) goto 10
	    call M$MALLOC (L_header_size,i)
	    call M$MOVE (hcb,%val(i),L_header_size)
	    Gx.nball = Gx.nball+1
	    Gx.lhsv(Gx.nball) = i
	    t1 = hcb.timecode + hcb.tstart
	    t2 = t1 + hcb.tdelta*(hcb.size-1)
	    if (Gx.nball .eq. 1) then
	      Gx.tmin = t1
	      Gx.tmax = t2
	    elseif (t1.gt.Gx.tmax .or. t2.le.Gx.tmin) then
	      call M$WARNING('State vector files do not overlap')
	    elseif (hcb.size .gt. 1) then
	      Gx.tmin = max (Gx.tmin, t1)
	      Gx.tmax = min (Gx.tmax, t2)
	    endif
  10	    str = str(j+1:)
	  enddo
	  Gx.time = max(Gx.tmin, min(Gx.tmax, Gx.time))
	  Gx.timel = -1

	elseif (type.eq.P_LAYER) then			! layer files
	  str = string
	  do while (M$LENGTH(str).gt.0)
	    j = M$SEARCH(str,'|')
	    if (j.eq.0) j = M$LENGTH(str)+1
	    hcb.file_name = str(1:j-1)
	    i = M$APICK_POST (hcb.file_name,j-1)
	    call M$MALLOC (L_layersz,ply)
	    call XGEO_GET_DATA (ly,hcb.file_name,0)
	    if (ly.nelem .ge. 0) then
	      Gx.nlayer = Gx.nlayer+1
	      Gx.lloc(Gx.nlayer) = ply
	    else
	      call MX$MESSAGE ('Couldnt open file: '//str(1:j-1),1.0)
	      call M$FREE (ply)
	    endif
	    str = str(j+1:)
	  enddo

	elseif (type.eq.P_EVENT) then			! event files
	  str = string
	  do while (M$LENGTH(str).gt.0)
	    j = M$SEARCH(str,'|')
	    if (j.eq.0) j = M$LENGTH(str)+1
	    hcb.file_name = str(1:j-1)
	    i = M$APICK_POST (hcb.file_name,j-1)
	    call XGEO_EVENTS_LOAD (hcb.file_name,0)
	    str = str(j+1:)
	  enddo

	elseif (type.eq.P_ALTDB) then			! altitude database
	  call M$INIT (hcb,string,'2000',' ',0)
	  call M$OPEN (hcb,HCBF_INPUT+HCBF_NOABORT)
	  if (.not. hcb.open) return
	  call M$FORCE1000 (hcb)
	  hcb.buf_type = 'D'
	  call M$MALLOC (L_header_size,Gx.lhadb)
	  call M$MOVE (hcb,%val(Gx.lhadb),L_header_size)
	endif

	return
	end


	subroutine XGEO_UNLOAD (type)
	include 'xgeoplot.inc'
	pointer (ply,ly)
	integer*4 type

	if (type.eq.P_SITE .or. type.eq.0) then	! site
	endif

	if (type.eq.P_SV .or. type.eq.0) then	! state vectors
	  do n = 1,Gx.nball
	    if (Gx.lhsv(n).ne.0) call M$CLOSE (%val(Gx.lhsv(n)))
	    call M$FREE(Gx.lhsv(n))
	  enddo
	endif

	if (type.eq.P_LAYER .or. type.eq.0) then	! layers
	  do n = Gx.nlayer,1,-1
	    ply = Gx.lloc(n) 
	    if (BTEST(ly.enable,PM_BS)) Gx.bsmask = 0
	    if (ly.lhcb.ne.0) call M$CLOSE( %val(ly.lhcb) )
	    call M$FREE (ly.bufloc)
	    call M$FREE (ly.lhcb)
	    call M$FREE (ply)
	  enddo
	  Gx.nlayer = 0	
	endif

	if (type.eq.P_EVENT .or. type.eq.0) then	! events
	  do n = 1,Gx.nevent
	    call M$FREE( Gx.eloc(n) )
	  enddo
	  Gx.nevent = 0
	endif

	if (type.eq.P_ALTDB .or. type.eq.0) then	! altitude dBase
	  if (Gx.lhadb.ne.0) call M$CLOSE (%val(Gx.lhadb))
	endif

	return
	end


	subroutine XGEO_RELOAD (type,index)
	include 'xgeoplot.inc'
	pointer (ply,ly)
	integer*4 type,index

	if (type.eq.P_SITE .or. type.eq.0) then	! site
	endif

	if (type.eq.P_SV .or. type.eq.0) then	! state vectors
	endif

	if (type.eq.P_LAYER .or. type.eq.0) then	! layers
	  do n = Gx.nlayer,1,-1
	    if (n.eq.index .or. index.eq.0) then
	      ply = Gx.lloc(n) 
	      if (ly.lhcb .eq. 0) then	! if not pipe
		if (BTEST(ly.enable,PM_BS)) Gx.bsmask = 0
		call M$FREE (ly.bufloc)
		call XGEO_GET_DATA (ly,ly.name,1)
	      endif
	    endif
	  enddo
	endif

	if (type.eq.P_ALTDB .or. type.eq.0) then	! altitude dBase
	endif

	return
	end


	subroutine XGEO_LIST (type,string,j,mode)
	include 'xgeoplot.inc'
	character*(*)  string
	record /HEADER/ hcb
	pointer (ply,ly)
	pointer (pev,ev)
	pointer (phcb,hcb)
	character seperator
	integer*4 type,mode,defenable

	if (mode.eq.1) seperator = ','
	if (mode.eq.2) seperator = '|'

	if (type.eq.P_SITE) then	! site
	endif

	if (type.eq.P_SV) then		! vehicles
	  j = 0
	  do n = 1,Gx.nball
	    phcb = Gx.lhsv(n)
	    call TRIM_NAME (hcb.file_name, string(j+1:) )
	    j = M$LENGTH(string)+1
	    string(j:j) = seperator
	  enddo
	endif

	if (type.eq.P_LAYER) then	! layers
	  j = 0
	  do k = mode-1,Gx.nlayer
	    ply = Gx.lloc(k)
	    if (mode.eq.1) then
	      string(j+1:) = '-- '//ly.name
	      if (BTEST(ly.enable,PM_GLOBAL)) then
		if (BTEST(ly.enable,PM_BASE)) string(j+1:j+1) = '+'
		if (BTEST(ly.enable,PM_ZOOM)) string(j+2:j+2) = '+'
	      endif
	    else
	      string(j+1:) = ly.name
	      defenable = P_ENABLE
	      if (ly.tunits.eq.1 .or. ly.type.eq.5000)
     &		defenable = IBCLR(defenable,PM_BS)
	      if (ly.enable .ne. defenable) then
		j = M$LENGTH(string)
		write(string(j+1:j+8),'(4h(EM=,i3.3,1h))') ly.enable
	      endif
	    endif
	    j = M$LENGTH(string)+1
	    string(j:j) = seperator
	  enddo
	endif

	if (type.eq.P_EVENT) then	! events
	  j = 0
	  do k = 1,Gx.nevent
	    pev = Gx.eloc(k)
	    if (mode.eq.1) then
	      string(j+1:) = '-- '//ev.name
	      if (BTEST(ev.enable,PM_GLOBAL)) then
		if (BTEST(ev.enable,PM_BASE)) string(j+1:j+1) = '+'
		if (BTEST(ev.enable,PM_ZOOM)) string(j+2:j+2) = '+'
	      endif
	    else
	      string(j+1:) = ev.name
	    endif
	    j = M$LENGTH(string)+1
	    string(j:j) = seperator
	  enddo
	endif

	return
	end



	subroutine XGEO_PLOT ()
	include 'xgeoplot.inc'
	call XGEO_RESET (0)
	call XGEO_REFRESH (0)
	return
	end

	subroutine XGEO_PIPES ()
	include 'xgeoplot.inc'
	pointer (ply,ly)
	jj = Gx.status
	do n = 1,Gx.nlayer
	  ply = Gx.lloc(n)
	  if (ly.pipemode.gt.PP_OFF .and. ly.piperead.ne.0) then
	    ii = Gx.status
	    call XGEO_GET_DATA (ly,' ',2)
	    if (Gx.status.ne.ii .and. BTEST(ly.enable,PM_BS)) 
     &			Gx.bsmask=JIBCLR(Gx.bsmask,n)
	  endif
	enddo
	if (Gx.status.ne.jj) call XGEO_REFRESH (0)
	return
	end


	subroutine XGEO_RESET (flag)
	include 'xgeoplot.inc'
	! flag: 0=normal -1=init 1=base 2=wide -2=narrow 3=resize
	real*4 ratio
c
c initialization pass
	if (flag .eq. -1) then	! Base reset
	  Gx.x1 = -180 + Gx.xoffs
	  Gx.x2 = +180 + Gx.xoffs
	  Gx.y1 = -90
	  Gx.y2 = +90
	  Gx.x1p = -1.2
	  Gx.x2p = +1.2
	  Gx.y1p = -1.2
	  Gx.y2p = +1.2
	  Mx.level = 1
	endif
	Gx.bsmask = 0
	Gx.timel = -1
c
c set text size
	m = min(8,Mx.width/102)
	call MX$SETFONT (m)
	if (flag.eq.3) return
c
c initialize plot request range
	k = Mx.level 
	Mx.level = 1
	if (flag .eq. 2 .or. flag.eq.-2) then	! wide/narrow view
	  dval = 6
	  if (flag.eq.2) dval = 1/dval
	  Mx.stk(1).xmin = Mx.stk(1).xmin * dval
	  Mx.stk(1).xmax = Mx.stk(1).xmax * dval
	  Mx.stk(1).ymin = Mx.stk(1).ymin * dval
	  Mx.stk(1).ymax = Mx.stk(1).ymax * dval
	elseif (Gx.view .eq. 1) then
	  if (flag.eq.1) then
	    Gx.x1 = Mx.stk(k).xmin
	    Gx.x2 = Mx.stk(k).xmax
	    Gx.y1 = Mx.stk(k).ymin
	    Gx.y2 = Mx.stk(k).ymax
	  endif
	  Mx.stk(1).xmin = Gx.x1
	  Mx.stk(1).xmax = Gx.x2
	  Mx.stk(1).ymin = Gx.y1
	  Mx.stk(1).ymax = Gx.y2
	else
	  if (flag.eq.1) then
	    Gx.x1p = Mx.stk(k).xmin
	    Gx.x2p = Mx.stk(k).xmax
	    Gx.y1p = Mx.stk(k).ymin
	    Gx.y2p = Mx.stk(k).ymax
	  endif
	  Mx.stk(1).xmin = Gx.x1p
	  Mx.stk(1).xmax = Gx.x2p
	  Mx.stk(1).ymin = Gx.y1p
	  Mx.stk(1).ymax = Gx.y2p
	  if (.not. Gx.prop .or. Gx.split) then
	    ! no adjustments
	  elseif (Mx.width .gt. Mx.height) then
	    ratio = float(Mx.width)/float(Mx.height)
	    Mx.stk(k).xmin = Mx.stk(k).xmin * ratio
	    Mx.stk(k).xmax = Mx.stk(k).xmax * ratio
	  else
	    ratio = float(Mx.height)/float(Mx.width)
	    Mx.stk(k).ymin = Mx.stk(k).ymin * ratio
	    Mx.stk(k).ymax = Mx.stk(k).ymax * ratio
	  endif
	endif
	return
	end



	subroutine XGEO_REFRESH (flag)
	! flag: 0=normal -1=nullbs 1=zoom 2=fast 3=fast/NoBS 4=justgrid
	include 'xgeoplot.inc'
	pointer (pev,ev)
c
c change perspective
	if (Gx.timel.ne.Gx.time .and. Gx.view.gt.2) then
	  call XGEO_SET_PERSP (Gx.view-2)	
	  Gx.bsmask = 0
	endif
	Gx.timel = Gx.time
c
c get new state vectors 
	do n = 1,Gx.nevent
	  pev = Gx.eloc(n)
	  if (ev.svmode .eq. PG_CURRENT) then
	    if (Gx.time .ne. ev.time .or. flag.le.0) 
     &		call XGEO_GET_SVS (ev, Gx.time)
	  elseif (ev.svmode .eq. PG_CURROFF) then
	    if (Gx.time+ev.dtime .ne. ev.time .or. flag.le.0) 
     &		call XGEO_GET_SVS (ev, Gx.time+ev.dtime)
	  endif
	enddo
c
c invalidate BS masks at all levels ?
	k = Mx.level
	if (Gx.bsmask.eq.0 .or. (Gx.zs.neqv.Gx.split) .or. flag.eq.-1) then
	  do i = 1,k
	    Gx.bsmap(i).mask = 0
	  enddo
	endif
	if (Gx.hold) return ! handle multiple messages before refresh
	if (Mx.icon) return ! dont draw when iconified
c
c split mode clears
	if (Gx.zs .and. Gx.zk.gt.1) call XGEO_DRAW_SZZ ()
	if (Gx.zs .neqv. Gx.split) call MX$CLEAR_WINDOW ()
	if (Gx.zk.gt.1 .and. k.eq.1) call MX$CLEAR_WINDOW ()
c
c draw view
	if (Gx.split .and. k.gt.1 .and. 
     &		(flag.le.0 .or. Gx.bsmask.eq.0) ) then
	  Mx.level = 1
	  call XGEO_REFRESH_VIEW (flag)
	  call XGEO_FACILITIES (2)
	  Mx.level = k
	endif
	call XGEO_REFRESH_VIEW (flag)
	call XGEO_FACILITIES (1)
c
c draw split zoom region rubber box
	if (Gx.split .and. k.gt.1) then
	  Gx.zx = (Mx.stk(k).xmin-Mx.stk(1).xmin)/Mx.stk(1).xscl+Mx.stk(1).x1
	  Gx.zy = (Mx.stk(1).ymax-Mx.stk(k).ymax)/Mx.stk(1).yscl+Mx.stk(1).y1
	  Gx.zw = (Mx.stk(k).xmax-Mx.stk(k).xmin)/Mx.stk(1).xscl
	  Gx.zh = (Mx.stk(k).ymax-Mx.stk(k).ymin)/Mx.stk(1).yscl
	  call XGEO_DRAW_SZZ ()
	endif
	Gx.zk = k
	Gx.zs = Gx.split
	Gx.bsmask = -1
c
c draw menu bar
	if (Gx.bar) then
	  call MX$TEXT (Mx.text_w,Mx.text_h,
     &		'Refresh   Viewpoint Functions Utilities Layers    '// 
     &		'Events    Options   KeyInfo')
	endif
c
c draw specs
	if (Gx.spec) call XGEO_SPECS ()
	Mx.stk(k+1) = Mx.stk(k)
	return
	end


	subroutine XGEO_REFRESH_VIEW (flag)
	include 'xgeoplot.inc'
	pointer (ply,ly),(pev,ev)
	record /LAYER/ lym
	record /T5REC/ p
	integer*4 curmask,newmask,addmask,delmask,bsmask
	integer*4 lowmask,drawmask
c
c set window borders
	k = Mx.level
	if (Gx.split .and. k.eq.1) then
	  ix1 = 0
	  ix2 = Mx.width/2
	elseif (Gx.split) then
	  ix1 = Mx.width/2
	  ix2 = Mx.width
	else
	  ix1 = 0
	  ix2 = Mx.width
	endif
	iy1 = 0
	iy2 = Mx.height
	if (Gx.bar) iy1 = iy1 + Mx.text_h
	if (Gx.spec) iy2 = iy2 - Mx.text_h*2

	if (Gx.axes) then
	  Mx.stk(k).x1 = ix1 + Mx.text_w*6
	  Mx.stk(k).x2 = ix2 - Mx.text_w*3
	  Mx.stk(k).y1 = iy1 + Mx.text_h*2
	  Mx.stk(k).y2 = iy2 - Mx.text_h*2
	else
	  Mx.stk(k).x1 = ix1 + 2
	  Mx.stk(k).x2 = ix2 - 4
	  Mx.stk(k).y1 = iy1 + 2
	  Mx.stk(k).y2 = iy2 - 4
	endif
	if (Gx.pan .and. k.gt.1 .and. .not.Gx.axes) then
	  Mx.stk(k).y1 = Mx.stk(k).y1 + Mx.text_h
	  Mx.stk(k).x2 = Mx.stk(k).x2 - Mx.text_h
	endif

	ix =  max(1,Mx.stk(k).x2-Mx.stk(k).x1)
	iy =  max(1,Mx.stk(k).y2-Mx.stk(k).y1)
	Mx.stk(k).xscl = (Mx.stk(k).xmax-Mx.stk(k).xmin)/ix
	Mx.stk(k).yscl = (Mx.stk(k).ymax-Mx.stk(k).ymin)/iy
	if (flag.eq.-9) return	! calculations for ZFILL done

	if (Gx.prop) then
	  dval =  Mx.stk(k).yscl / Mx.stk(k).xscl 
	  if (dval .gt. 1.d0) then
	    Mx.stk(k).xscl = Mx.stk(k).yscl
	    if (k.eq.1 .and. Gx.view.gt.1) then
	      dval = (Mx.stk(k).xmax+Mx.stk(k).xmin)*.5
	      Mx.stk(k).xmin = dval-ix*Mx.stk(k).xscl*.5
	      Mx.stk(k).xmax = dval+ix*Mx.stk(k).xscl*.5
	    else
	      i = (Mx.stk(k).x2-Mx.stk(k).x1)/2
	      ii = (Mx.stk(k).x2+Mx.stk(k).x1)/2
	      Mx.stk(k).x1 = ii - i/dval
	      Mx.stk(k).x2 = ii + i/dval
	    endif
	  else
	    Mx.stk(k).yscl = Mx.stk(k).xscl
	    if (k.eq.1 .and. Gx.view.gt.1) then
	      dval = (Mx.stk(k).ymax+Mx.stk(k).ymin)*.5
	      Mx.stk(k).ymin = dval-iy*Mx.stk(k).yscl*.5
	      Mx.stk(k).ymax = dval+iy*Mx.stk(k).yscl*.5
	    else
	      i = (Mx.stk(k).y2-Mx.stk(k).y1)/2
	      ii = (Mx.stk(k).y2+Mx.stk(k).y1)/2
	      Mx.stk(k).y1 = ii - i*dval
	      Mx.stk(k).y2 = ii + i*dval
	    endif
	  endif
	endif
c
c initialize transformation variables
	call DRAW_INIT()
c
c invoke clip region
	call MX$CLIP (Mx.stk(k).x1, Mx.stk(k).y1,
     &    Mx.stk(k).x2-Mx.stk(k).x1, Mx.stk(k).y2-Mx.stk(k).y1)
c
c backing store layer masks
	curmask = Gx.bsmap(k).mask
	if (Gx.bsmap(k).ip.eq.0 .or. flag.ne.0) curmask = 0
	lowmask = 0
	newmask = 0
	bsmask = 0
	do n = 0,Gx.nlayer
	  ply = Gx.lloc(n)
	  if (XGEO_ENABLED(ly.enable,flag) .and.
     &		ly.nelem.gt.0) newmask = JIBSET (newmask,n)
	  if (BTEST(ly.enable,PM_BS)) bsmask = JIBSET (bsmask,n)
	  if (ly.type .eq. 2000) lowmask = JIBSET (lowmask,n)
	enddo
	if (flag.eq.4) newmask = 1	! grid layer only
c--CODE CHANGE: (.not. to not) for Intel Compiler JLR/EHG
	delmask = curmask .and. not(newmask) 
	if (delmask .ne. 0) curmask = 0
	addmask = newmask .and. not(curmask) 
c--	
	if (curmask.ne.0) then
	  call MX$MAP (Gx.bsmap(k),0,0,0,0,L_PutMap)
	  call MX$CLEAR_AREA (ix1,iy1,ix2-ix1,Mx.stk(k).y1-iy1)
	  call MX$CLEAR_AREA (ix1,Mx.stk(k).y2,ix2-ix1,iy2-Mx.stk(k).y2)
	  call MX$CLEAR_AREA (ix1,iy1,Mx.stk(k).x1-ix1,iy2-iy1)
	  call MX$CLEAR_AREA (Mx.stk(k).x2,iy1,ix2-Mx.stk(k).x2,iy2-iy1)
	else
	  call MX$CLEAR_AREA (ix1,iy1,ix2-ix1,iy2-iy1)
	endif
c
c draw pan bars
	ii = XGEO_PAN_HANDLE (0,i) 
c
c plot layers
	drawmask = addmask .and. lowmask .and. bsmask
	do n = 0,Gx.nlayer
	  if (BTEST(drawmask,n)) call XGEO_DRAW_LAYER (%val(Gx.lloc(n)))
	enddo
c--CODE CHANGE: (.not. to not) for Intel Compiler JLR/EHG	
	drawmask = addmask .and. not(lowmask) .and.  bsmask 
c--	
	do n = 0,Gx.nlayer
	  if (BTEST(drawmask,n)) call XGEO_DRAW_LAYER (%val(Gx.lloc(n)))
	enddo
	drawmask = addmask .and. bsmask
	if (Gx.bs .and. drawmask.ne.0 .and. flag.le.2) ! Update BS pixmap
     +	    call MX$MAP (Gx.bsmap(k),Mx.stk(k).x1,Mx.stk(k).y1,
     +		Mx.stk(k).x2-Mx.stk(k).x1,Mx.stk(k).y2-Mx.stk(k).y1,L_NewMap)
c .. add line thickness control (from THICK==n keyword) for non-BS layers,
c .. which is typically how user-generated map overlay files come in
c--CODE CHANGE: (.not. to not) for Intel Compiler JLR/EHG
	drawmask = addmask .and. not(bsmask) 
c--	
	do n = 0,Gx.nlayer
	  ply = Gx.lloc(n) 
	  if (BTEST(drawmask,n)) then
	    call MX$DRAW_MODE (ly.thick)
	    call XGEO_DRAW_LAYER (ly)
	    call MX$DRAW_MODE (1)
	  endif
	enddo
	if (flag.gt.2 .or. .not. Gx.bs) newmask = 0
	Gx.bsmap(k).mask = newmask .and. bsmask
c
c plot events
	if (Gx.evnt) then
	  do n = 1,Gx.nevent
	    pev = Gx.eloc(n)
	    if (XGEO_ENABLED(ev.enable,flag)) call XGEO_DRAW_EVENT (n)
	  enddo
	  if (Gx.cep) call XGEO_DRAW_EVENT (0)
	endif
c
c plot vehicle cones
	if (Gx.vcon) call XGEO_VEHICLE_CONE (0)
c
c plot marker
	if (Gx.mark .and. 
     &	   (Gx.mrk.x.ne.0 .or. Gx.mrk.y.ne.0 .or. Gx.mrk.z.ne.0)) then
	  lym.type = 5000
	  lym.name = '*'
	  lym.nelem = 1
	  lym.color = -1
	  lym.bufloc = %loc(Gx.mrk)
	  lym.enable = P_ENABLE
	  call XGEO_DRAW_LAYER (lym)
	endif
c
c clear clip
	call MX$CLIP (0,0,0,0)
c
c draw axes
	if (Gx.view .eq. 1) then
	  if (Gx.axes) call MX$DRAWAXIS (5, 5, 61, 60, 0)
	else
	  if (Gx.axes) call MX$DRAWAXIS (5, 5, 0, 0, 0)
	endif
c
	return
	end


	subroutine XGEO_DRAW_SZZ ()
	include 'xgeoplot.inc'
	call MX$DRAW_BOX (-2,Gx.zx,Gx.zy,Gx.zw,Gx.zh)
	i = 10
	ix = Gx.zx + Gx.zw/2
	iy = Gx.zy + Gx.zh/2
	call MX$DRAW_LINE (-2,ix,iy-i,ix,iy+i)
	call MX$DRAW_LINE (-2,ix-i,iy,ix+i,iy)
	return
	end



	subroutine XGEO_EVENT ()
	include 'xgeoplot.inc'
	pointer (pev,ev),(ply,ly)
	external XGEO_SPECS, XGEO_MAINMENU
	real*4 fdif

	if (Gx.rtdelta .gt. 0.0) then
	  fval = SECNDS(0.0)
	  fdif = fval - Gx.rtlast		! <0 if new day
	  if (fdif .ge. Gx.rtdelta .or. fdif .lt. 0) then
	    Gx.rtlast = INT(fval/Gx.rtdelta)*DBLE(Gx.rtdelta)
	    call M$NOW (Gx.time, dval)
	    call M$DRSLT ('G_TIME',Gx.time)
	    call XGEO_GET_SVS (%val(Gx.eloc(1)),Gx.time)
	    call XGEO_REFRESH (0)
	  endif
	endif

	do while ( MX$IFEVENT() )
	 Gx.status = Gx.status + 1

	if ( XGEO_PAN_HANDLE(1,i) ) then
	  if (i.gt.0) then
	    call XGEO_REFRESH (2)
	    call M$PAUSE (0.2)
	    ii = MX$IFEVENT()
	  endif
	  m = 3
	  do while (Mx.button_release.eq.0)
	    ii = XGEO_PAN_HANDLE(2,i)
	    if (i.gt.0) then
	      fval = SECNDS(0.0)
	      call XGEO_REFRESH (m)
	      ii = MX$IFEVENT()
	      fval = SECNDS(fval)
	      if (fval.lt.Gx.speed) call M$PAUSE (Gx.speed-fval)
	      if (fval.gt.0.5) m = 4
	    else
	      call MX$EVENT()
	    endif
	  enddo
	  call XGEO_REFRESH (0)
	  ii = XGEO_PAN_HANDLE(3,i)

	elseif (Mx.motion_notify) then
	  call XGEO_SPECS ()
	  if (Gx.cntrls.eq.2) then
	    call M$DRSLT (Mx.xtag,Gx.ret.x)
	    call M$DRSLT (Mx.ytag,Gx.ret.y)
	    call M$DRSLT (Mx.ztag,Gx.ret.z)
	  endif

	elseif (Mx.button_press.ne.0 .and.		! Time Functions
     &	  Gx.spec .and. Mx.ypos .gt. Mx.height-2*Mx.text_h .and. 
     &    Mx.xpos .gt. Mx.width-15*Mx.text_w) then

	  ! Alter time fuction replaces animate mode
	  if (Mx.button_press .eq. 2) then
	    call XGEO_GET_TIME ('TIME = ',Gx.time)
	    call XGEO_GET_SVS (%val(Gx.eloc(1)),Gx.time)
	    call XGEO_REFRESH (0)
	  else
	    i = (Mx.width-Mx.xpos)/Mx.text_w
	    dval = 3600
	    if (i .le. 6) dval = 600
	    if (i .le. 4) dval = 60
	    if (i .le. 2) dval = 1
	    if (Mx.ypos .lt. Mx.height-Mx.text_h) dval = 86400
	    if (Mx.button_press .eq. 3) dval = -dval
  22	    if (dval.lt.0.d0 .and. Gx.time.le.Gx.tmin) goto 23
	    if (dval.gt.0.d0 .and. Gx.time.ge.Gx.tmax) goto 23
	    Gx.time = max(Gx.tmin, min(Gx.tmax, Gx.time+dval))
	    call XGEO_GET_SVS (%val(Gx.eloc(1)),Gx.time)
	    call XGEO_REFRESH (0)
	    if (Mx.button_press.ne.0) then
	      i = MX$IFEVENT()
	      if (Mx.button_release.ne.0) goto 23
	      call M$PAUSE (0.2)
	    endif
	    call M$PAUSE (Gx.speed)	! max animate speed
	    do while (MX$IFEVENT() )
	      if (Mx.button_release.ne.0) goto 23
	    enddo
	    goto 22
  23	    continue
	  endif
	  call M$DRSLT ('G_TIME',Gx.time)

	elseif (Mx.button_press.eq.2 .or. Mx.key_upchar.eq.'M') then

	  m = MX$MENUPLUS (0,'X-GEOPLOT;Refresh,Viewpoint,Functions,'//
     &		'Layers,Vehicles,Events,Conics,Utilities,Options,'//
     &		'Key Info,Exit', 0, L_CBenable, XGEO_MAINMENU)

	elseif (Mx.button_press.eq.1) then

	  if (Gx.bar .and. Mx.xpos.lt.Mx.text_w*76 .and. 
     &		Mx.ypos.lt.Mx.text_h*3/2) then			! menu bar
	    m = max(1,min(8,(Mx.xpos/Mx.text_w+9)/10))
	    call MAINMENU (m)
	  endif

	  call MX$RUBBERBOX (ix1,iy1,ix2,iy2,XGEO_SPECS)
	  if (ix2-ix1.lt.5 .and. iy2-iy1.lt.5) then		! set origin
	    call XGEO_SPECS ()
	    if (Gx.sorg) then
	      Gx.mrk.x = Gx.ret.x
	      Gx.mrk.y = Gx.ret.y
	      Gx.mrk.z = Gx.ret.z
	      call MSGDO (MSK_MARK, Mx.level)
	    endif
	    if (Gx.cntrls.eq.1) then
	      call M$DRSLT (Mx.xtag,Gx.ret.x)
	      call M$DRSLT (Mx.ytag,Gx.ret.y)
	      call M$DRSLT (Mx.ztag,Gx.ret.z)
	    endif
	    call XGEO_REFRESH (0)
	  elseif (ix2-ix1.gt.5 .and. iy2-iy1.gt.5) then		! zoom
	    call XGEO_ZOOM (ix1,iy1,ix2,iy2)
	  endif

	elseif (Mx.button_press.eq.3) then	
	  k = Mx.level
	  if (k.gt.1) then			! unzoom
	    if (Gx.bsmap(k).ip.ne.0) 
     &		call MX$MAP (Gx.bsmap(k),0,0,0,0,L_FreeMap)
	    Mx.level = k-1
	    call XGEO_REFRESH (0)
	    call MSGDO (MSK_UNZOOM, Mx.level)
	  endif

	elseif (Mx.key_upchar.eq.'A' .and. Gx.nball.gt.0) then! toggle abscissa
	  Gx.iballa = mod(Gx.iballa,Gx.nball)+1
	  Gx.iballb = mod(Gx.iballb,Gx.nball)+1
	  if (Gx.nevent.gt.0) then
	    pev = Gx.eloc(1)
	    ev.iballa = Gx.iballa
	    ev.iballb = Gx.iballb
	    call XGEO_GET_SVS (ev,Gx.time)
	    call XGEO_REFRESH (0)	    
	  endif

	elseif (Mx.key_upchar .eq. 'G') then		! geolocate
	  if (XGEO_LOCATE().gt.0) then
	    call M$DRSLT ('G_LON',Gx.mrk.lon)
	    call M$DRSLT ('G_LAT',Gx.mrk.lat)
	    call M$DRSLT ('G_ALT',Gx.mrk.alt)
	    call XGEO_REFRESH (0)	    
	  else
	    call MX$MESSAGE ('Could not converge',-1.0)
	  endif

	elseif (Mx.key_upchar .eq. 'O') then		! set origin
	  call XGEO_SPECS ()
	  Gx.mrk.x = Gx.ret.x
	  Gx.mrk.y = Gx.ret.y
	  Gx.mrk.z = Gx.ret.z
	  call XGEO_REFRESH (0)

	elseif (Mx.key_upchar.eq.'B') then		! boxes facility
	  call XGEO_FACILITIES (-3)
	elseif (Mx.key_upchar.eq.'P') then		! points facility
	  call XGEO_FACILITIES (-4)

	elseif (Mx.key_char.eq.'w') then		! wide view 
	  call XGEO_RESET (-2)
	  call XGEO_REFRESH (0)
	elseif (Mx.key_char.eq.'W') then		! narrow view 
	  call XGEO_RESET (2)
	  call XGEO_REFRESH (0)

	elseif (Mx.key_upchar.eq.'U') then		! toggle units
	  Gx.units = mod(Gx.units,3) + 1
	  if (Gx.spec) call XGEO_SPECS ()

	elseif (Mx.key_upchar.eq.'V') then		! toggle view 
	  if (Mx.key_char.eq.'V') then 
	    Gx.view = Gx.view - 1
	    if (Gx.view .le. 0) Gx.view = Gx.nball+2
	  else
	    Gx.view = Gx.view + 1
	    if (Gx.view .gt. Gx.nball+2) Gx.view = 1
	  endif
	  call XGEO_RESET (0) 
	  call XGEO_REFRESH (0) 

	elseif (Mx.key_upchar .eq. 'E') then		! toggle events
	  Gx.evnt = .not. Gx.evnt
	  call XGEO_REFRESH (0)

	elseif (Mx.key_upchar .eq. 'L') then		! layers hot-key
	  call MAINMENU (4)

	elseif (Mx.key_press .ge. 48+1 .and. 
     &		Mx.key_press .le. 48+Gx.nevent ) then
	  k = Mx.key_press - 48 
	  call XGEO_EVENT_DIALOG (k)

	elseif (Mx.key_upchar .eq. 'H') then		! hardcopy
	  str = 'xgeoplot'
	  if (MX$HARDCOPY_OPEN(str, -1)) then
	    call XGEO_REFRESH (-1)
	    call MX$HARDCOPY_CLOSE()
	  endif

	elseif (Mx.key_press.eq.13 .and. Mx.noexit.le.0) then   ! the END
	  Gx.status = -1
	  return
	endif

	if (Mx.resize) then
	  call XGEO_RESET (3)
	  call XGEO_REFRESH (0)
	elseif (Mx.exposure) then
	  call XGEO_REFRESH (0)
	endif

	enddo
	return
	end

	function XGEO_ENABLED (enables,flag)
	include 'xgeoplot.inc'
	logical*4 ok
	integer*4 enables
	ok = BTEST(enables,PM_GLOBAL)
	if (flag.ge.2) then
	  if (Mx.level.eq.1 .and. .not.BTEST(enables,PM_FBASE)) ok=0
	  if (Mx.level.gt.1 .and. .not.BTEST(enables,PM_FZOOM)) ok=0
	else
	  if (Mx.level.eq.1 .and. .not.BTEST(enables,PM_BASE) ) ok=0
	  if (Mx.level.gt.1 .and. .not.BTEST(enables,PM_ZOOM) ) ok=0
	endif
	if (flag .eq. -9) then
	  enables = MX$MENUPLUS (enables,'ENABLES;Global,Labels,'//
     &	  'BackStore,Base,Zoom,FastBase,FastZoom', -1, 0, 0)
	endif
	XGEO_ENABLED = ok
	return
	end

	subroutine MAINMENU (item)
	include 'xgeoplot.inc'
	record /CBSTRUCT/ cb
	integer*4 item
	cb.item = item
	call XGEO_MAINMENU (cb)
	return
	end


	subroutine XGEO_MAINMENU (cb)
	include 'xgeoplot.inc'
	record /CBSTRUCT/ cb
	pointer (pev,ev),(ply,ly)
	external XGEO_SPECS, XGEO_SETTINGS
	integer*4 optmask

	m = cb.item
	cb.status = L_CBokaction	! OK status

	  if (m.eq.1) then				! refresh
	    call XGEO_REFRESH (-1)
	    cb.status = L_CBexit

	  elseif (m.eq.2) then				! viewpoint
	    n = Gx.view
	    Gx.view = MX$MENU (min(n,5),'VIEWPOINT;'//
     &		'Lat/Lon,Perspective,Vehicle 1,Vehicle 2,Vehicle N')
	    dval = Gx.pdis
	    if (Gx.view.eq.2) call MX$DPROMPT 
     &		('Perspective distance (Earth Radii): ',Gx.pdis)
	    Gx.pdis = max(1.01d0,Gx.pdis)
	    if (Gx.view.eq.5) then
	      i = max(1,n-2)
	      if (Gx.nball.gt.0) call MX$LVAL ('N-Vehicle: ',i,1,Gx.nball,1,0)
	      Gx.view = i+2
	    endif 
	    if (Gx.view .gt. Gx.nball+2) then
	      call MX$MESSAGE ('Vehicle does not exist',-1.0)
	      Gx.view = Gx.nball+2
	    endif
	    if (Gx.view.ne.n .or. dval.ne.Gx.pdis) then
	      call XGEO_RESET (0)
	      call XGEO_REFRESH (0)
	    endif

	  elseif (m.eq.3) then				! Functions
	    n = MX$MENU (0,'FUNCTIONS;Event Region,Base Region,'//
     +		'Base Reset,Back Off,Select Locate Mode,Set STD,'//
     +		'Dump Predicts,Set Marker,Grid Size,Spin,RealTime')
	    if (n.eq.1) then				! Event Region
	      call MX$MESSAGE ('Zoom on event region',-1.0)
	      call MX$RUBBERBOX (ix1,iy1,ix2,iy2,XGEO_SPECS)
	      Mx.xpos = ix1
	      Mx.ypos = iy1
	      call XGEO_SPECS ()
	      Gx.rlon1 = Gx.ret.lon
	      Gx.rlat2 = Gx.ret.lat
	      Mx.xpos = ix2
	      Mx.ypos = iy2
	      call XGEO_SPECS ()
	      Gx.rlon2 = Gx.ret.lon
	      Gx.rlat1 = Gx.ret.lat
	      call XGEO_REFRESH (0)
	    elseif (n.eq.2) then			! Base Region
	      call XGEO_RESET (1)
	      call XGEO_REFRESH (0)
	    elseif (n.eq.3) then			! Base Reset
	      call XGEO_RESET (-1)
	      call XGEO_REFRESH (0)
	    elseif (n.eq.4) then			! Back Off
	      call XGEO_RESET (-2)
	      call XGEO_REFRESH (0)
	    elseif (n.eq.5) then			! Select Locate Mode
	      Gx.glmode = MX$MENU(Gx.glmode,
     +		'GEOLOCATION MODE;Tdoa/Fdoa,Tdoa/Tdoa,Fdoa/Fdoa')
	      j = 0
	      do k = 1,Gx.nevent
		pev = Gx.eloc(k)
		str(j+1:) = ev.name
		j = M$LENGTH(str)+1
		str(j:j) = ','
	      enddo
	      Gx.glevn1 = MX$MENU (Gx.glevn1,'Event #1;'//str(1:j-1))
	      Gx.glevn2 = MX$MENU (Gx.glevn2,'Event #2;'//str(1:j-1))
	    elseif (n.eq.6) then			! set STD
	      if (gx.glmode.eq.1) then
	        call MX$FPROMPT ('Tdoa-1 STD (sec): ',Gx.glstd1)
	        call MX$FPROMPT ('Fdoa-2 STD (Hz): ',Gx.glstd2)
	      elseif (gx.glmode.eq.2) then
	        call MX$FPROMPT ('Tdoa-1 STD (sec): ',Gx.glstd1)
	        call MX$FPROMPT ('Tdoa-2 STD (sec): ',Gx.glstd2)
	      elseif (gx.glmode.eq.3) then
	        call MX$FPROMPT ('Fdoa-1 STD (Hz): ',Gx.glstd1)
	        call MX$FPROMPT ('Fdoa-2 STD (Hz): ',Gx.glstd2)
	      endif
	      call MX$FPROMPT ('CEP Contour step (m): ',Gx.cepstep)
	      call XGEO_REFRESH (0)
	    elseif (n.eq.7) then			! write predicts 
	      i = MX$MENU (0,'PREDICTS;Single,Time Range')
	      dval = 300
	      if (i.eq.1) then
	        t1 = Gx.time
	        call XGEO_GET_TIME ('Start Time = ',t1)
		t2 = t1	
	      elseif (i.eq.2) then
	        t1 = Gx.tmin
	        t2 = Gx.tmax
	        call XGEO_GET_TIME ('Start Time = ',t1)
	        call XGEO_GET_TIME ('End Time = ',t2)
	        call XGEO_GET_TIME ('Delta Time = ',dval)
	      endif
	      if (i.ne.0) then
	        str = ' '
		call MX$PROMPT ('Filename: ',str(1:40))
	        call XGEO_DUMP_PREDICTS (%val(Gx.eloc(1)),
     +		  t1,t2,dval,str(1:40))
	      endif
	    elseif (n.eq.8) then			! set Marker
	      call MX$DPROMPT ('Marker Latitude: ',Gx.mrk.lat)
	      call MX$DPROMPT ('Marker Longitude: ',Gx.mrk.lon)
	      call MX$DPROMPT ('Marker Altitude: ',Gx.mrk.alt)
	      call MX$PROMPT  ('Marker Name: ',Gx.mrk.name)
	      Gx.mrk.name_len = M$LENGTH(Gx.mrk.name)
	      call MX$PROMPT  ('Marker Text: ',Gx.mrk.text)
	      Gx.mrk.text_len = M$LENGTH(Gx.mrk.text)
	      i = Gx.mrk.color
	      call MX$LPROMPT  ('Marker Color: ',i)
	      Gx.mrk.color = i
	      call XGEO_REFRESH (0)
	    elseif (n.eq.9) then			! set Grid resolution
	      call MX$LPROMPT ('Grid Resolution: ',Gx.nmg)
	      Gx.nmg = max(3,min(101,Gx.nmg))
	      call M$REMALLOC (Gx.nmg*Gx.nmg*4*5,Gx.lmg)
	      call XGEO_REFRESH (0)
	    elseif (n.eq.10) then			! Spin
	      n = 1
	      if (Gx.view .eq. 1) then
		Gx.view = 2 
		call XGEO_RESET (0)
	      endif
	      ii = 5
	      i = Gx.bs
	      Gx.bs = .false.
	      call MX$SETBAR ('Spin Mode')
	      do while (n.ne.3)
		if (n.eq.2) then
		  call M$PAUSE (Mc.pause)
		else
		  Gx.plon = Gx.plon + ii
		  if (Gx.plon .gt. 180.+Gx.xoffs) Gx.plon = Gx.plon - 360
		  call XGEO_REFRESH (0)
		endif
	        if ( MX$IFEVENT() ) then
		  if (Mx.button_press.eq.2) then
	            n = MX$MENU (n,'SPIN;Go,Pause,Stop,Delta')
		    if (n.eq.4) call MX$LPROMPT ('Longitude division: ',ii)
		  endif
		endif
	      enddo
	      Gx.bs = i
	      call MX$UNSETBAR ()
	    elseif (n.eq.11) then			! Real Time
	      call MX$FVAL ('RT Delta (sec): ',Gx.rtdelta, 
     +                       -1.0,300.0, 1.0, 0)
	    endif

	  elseif (m.eq.4) then				! layers
  17	    ix = Mx.xpos
	    iy = Mx.ypos
	    str = 'LAYERS;Newfile,'
	    j = M$LENGTH(str)
	    call XGEO_LIST (P_LAYER,str(j+1:),j,1)
	    j = M$LENGTH(str)
	    Mx.xpos = ix
	    Mx.ypos = iy
	    n = MX$MENU (0,str(1:j))
	    if (n.eq.1) then		! newfile
	      str = ' '
	      if (MX$FILES(str(1:80)) .gt. 0) then
		j = M$CLEAN(str,str)-1
		call XGEO_LOAD (P_LAYER,str(1:j))
	        Mx.exposure = .true.
	      endif
	    elseif (n.gt.1) then		! modify layer
	      k = n-2
	      call XGEO_LAYER_MOD (k)
	      Mx.exposure = .true.
	    endif	
	    if (n.ne.0) goto 17
	    if (Mx.exposure) call XGEO_REFRESH (0)

	  elseif (m.eq.5) then				! vehicles
  18	    ix = Mx.xpos
	    iy = Mx.ypos
	    str = 'VEHICLES;- Enable,Newfile,'
	    j = M$LENGTH(str)
	    call XGEO_LIST (P_SV,str(j+1:),j,1)
	    j = M$LENGTH(str)
	    if (Gx.vcon) str(10:10) = '+'
	    Mx.xpos = ix
	    Mx.ypos = iy
	    n = MX$MENU (0,str(1:j))
	    if (n.eq.1) then		! toggle
	      Gx.vcon = .not. Gx.vcon
	      call XGEO_REFRESH (0)
	    elseif (n.eq.2) then	! newfile
	      str = ' '
	      if (MX$FILES(str(1:80)) .gt. 0) then
		j = M$CLEAN(str,str)-1
		call XGEO_LOAD (P_SV,str(1:j))
	      endif
	    elseif (n.gt.2) then		! modify vehicle
	      k = n-2
	      i = MX$MENU(0,'VEHICLE;ReOpen,Delete')
	      if (i.eq.1) then
		call M$HCBFUNC (%val(Gx.lhsv(k)),HCBF_REOPEN,0)
	      elseif (i.eq.2) then
		call M$CLOSE (%val(Gx.lhsv(k)))
		do k = k,Gx.nball-1
		  Gx.lhsv(k) = Gx.lhsv(k+1)
		enddo
		Gx.nball = Gx.nball-1
	      endif
	    endif	
	    if (n.ne.0) goto 18

	  elseif (m.eq.6) then				! events
	    str = 'EVENTS;- Enable,Save Events,Load Events,New Event,'
	    j = M$LENGTH(str)
	    call XGEO_LIST (P_EVENT,str(j+1:),j,1)
	    j = M$LENGTH(str)
	    if (Gx.evnt) str(8:8) = '+'
	    n = MX$MENU (0,str(1:j))
	    if (n.eq.1) then				! Toggle Enable
	      Gx.evnt = .not. Gx.evnt
	      call XGEO_REFRESH (0)
	    elseif (n.eq.2) then			! Save event
	      str = ' '
	      call MX$PROMPT ('File name: ',str(1:40))
	      call XGEO_EVENTS_SAVE (str,0)
	    elseif (n.eq.3) then			! Load event
	      str = ' '
	      call MX$PROMPT ('File name: ',str(1:40))
	      call XGEO_EVENTS_LOAD (str,0)
	      call XGEO_REFRESH (0)
	    elseif (n.eq.4) then			! new event
	      call M$MALLOC (L_eventsz,pev)
	      if (Gx.nevent.eq.0) then
	        if (Gx.nball.eq.0) then
		  call MX$MESSAGE('No vehicles present',1.0)
		  return
		endif
	        Gx.nevent = Gx.nevent+1
	        Gx.eloc(Gx.nevent) = pev
	        call XGEO_EVENT_INIT (Gx.nevent)
	      else
	        call M$MOVE (%val(Gx.eloc(Gx.nevent)),ev,L_eventsz)
	        Gx.nevent = Gx.nevent+1
	        Gx.eloc(Gx.nevent) = pev
	        call MX$MESSAGE('Duplicated last event for modify',-1.0)
	      endif
	      call XGEO_EVENT_MOD (Gx.nevent)
	      call XGEO_REFRESH (0)
	    elseif (n.gt.4) then			! modify event
	      k = n-4
	      pev = Gx.eloc(k)
	      call XGEO_EVENT_MOD (k)
	      call XGEO_REFRESH (0)
	    endif	

	  elseif (m.eq.7) then				! conics
	    ix = Mx.xpos
	    iy = Mx.ypos
 27	    Mx.xpos = ix
	    Mx.ypos = iy
	    n = MX$MENU(0,'CONIC;Mode,Vehicle,Divisions,Levels,Angle')
	    if (n.eq.1) then
	      Gx.conemod = MX$MENU(Gx.conemod+1,
     &			'CONE MODE;Off,Slice,Centerline,Both') - 1
	    elseif (n.eq.2) then
	      Gx.coneveh = MX$MENU(Gx.coneveh+2,
     &			'CONE VEHICLE;All,None,1,2,3,4,5,6,7,8') - 2
	    elseif (n.eq.3) then
	      call MX$LVAL ('Cone Divisions: ',Gx.conediv, 2, 240, 1, 0)
	    elseif (n.eq.4) then
	      call MX$LVAL ('Cone Levels: ',Gx.conelev, -1, 16, 1, 0)
	    elseif (n.eq.5) then
	      call MX$FVAL ('Cone Angle: ',Gx.coneang, 
     &                       0.0, 90.0, 0.1, 0)
	    endif
	    call XGEO_REFRESH (0)
	    if (n.ne.0) goto 27

	  elseif (m.eq.8) then				! utilities
	    n = MX$MENU (0,'UTILS;'//
     &	'Colormap,Cntrls,Boxes,Points,Annotate,Speed,Units,Linewidth')
	    if (n.eq.1) then
	      call XGEO_COLORS (-1,' ')
	    elseif (n.eq.2) then
	      n = Gx.cntrls
	      Gx.cntrls = MX$MENU (Gx.cntrls+3, 'CONTROLS OPTIONS;'//
     &	'-Continuous,-LM Click, Off,+LM Click,+Continuous, Assign,') - 3
	      if (Gx.cntrls .eq. 3) then
		Gx.cntrls = n
		call MX$PROMPT ('Mouse Longitude -> ', Mx.xtag)
		call M$UPPERCASE (Mx.xtag)
		call MX$PROMPT ('Mouse Latitude -> ', Mx.ytag)
		call M$UPPERCASE (Mx.ytag)
		call MX$PROMPT ('Mouse Altitude -> ', Mx.ztag)
		call M$UPPERCASE (Mx.ztag)
	      endif
	    elseif (n.ge.3 .and. n.le.5) then
	      call XGEO_FACILITIES (n)
	    elseif (n.eq.6) then
	      j = NINT(1.0/Gx.speed)
	      call MX$LVAL ('Max frame rate (Hz): ',j,1,50,1,0)
	      Gx.speed = 1./float( max(1,j) )
	    elseif (n.eq.7) then
	      Gx.units = MX$MENU(Gx.units,
     &		'SPEC UNITS;Metric,English,Nautical')
	      if (Gx.spec) call XGEO_SPECS ()
	    elseif (n.eq.8) then
	      Gx.linew = MX$MENU(Gx.linew,
     &		'Event Line Width;1,2,3,4,5')
	      call XGEO_REFRESH (0)
	    endif

	  elseif (m.eq.9) then				! Options

	    optmask = 0
	    do i = 1,19
	      if (Gx.options(i)) optmask = IBSET(optmask,i)
	    enddo
	    optmask = MX$MENUPLUS (optmask,'OPTIONS;Labels,Axes,Specs,'//
     &		'Aspect,Panning,Split,Menu,BStore,Events,CEP,SetOrg,'//
     &		'Marker,VehCone,Zfill,Hold,SV=Hor,SView,MView,AoASpec',
     &		-1,L_CBenable,XGEO_SETTINGS)

	  elseif (m.eq.10) then					! keypress info
	    call MX$MHELP ('X-GEOPLOT','Key Info')

	  elseif (m.eq.11) then					! exit
	    if (Mx.noexit.gt.1) then
	      call MX$MESSAGE ('Exit disabled', 1.0)
	    else
	      Gx.status = -1
	      cb.status = L_CBexit
	    endif
	  endif

	return
	end


	subroutine XGEO_SETTINGS (cb)
	include 'xgeoplot.inc'
	record /CBSTRUCT/ cb
	Gx.options(cb.item) = .not. Gx.options(cb.item)
	if (cb.item.eq.6 .and. Mx.level.eq.1) call XGEO_RESET (0)
	if (cb.item.ge.2 .and. cb.item.le.7) then
	  call XGEO_REFRESH (-1)
	else
	  call XGEO_REFRESH (0)
	endif
	cb.status = L_CBokaction
	return
	end



	subroutine XGEO_MESSAGE ()
	include 'xgeoplot.inc'
	include 'qmessdef.inc'
	include 'qmessdecl.inc'
	pointer (pev,ev),(ply,ly)
	byte options(16)
	character*8 action
	integer*4 sid,color,msgid,mimic,info,lastmode
	record /T5REC/ t5
	pointer (pt5,t5)

	do while (M$GET_MSG(MQH,MQD,L_mqd,0.0,0) .eq. 1)
	  Gx.status = Gx.status + 1

	  Gx.xf.mqdloc.p = %loc(MQD)
	  Gx.xf.iarg = MQH.info
	  sid = MQH.sender
	  info = MQH.info
	  msgid = Gx.xf.msgid
	  mimic = Gx.mimic

	if (MX$GRAPHICS_MSG(Mx,MQH,MQD,msgid,Gx.xf.msgmask).ne.0) then
	  if (Mx.breakflag) then
	    Gx.status = -1
	    return
	  endif

	elseif (MX$MIMIC_MSG(MQH,MQD,mimic).lt.0) then
	  ! info message but not from who we're mimicing

	elseif (MQH.name .eq. '?ACK') then
	  call M$SEND_MSGL ('=ACK',0,sid,Gx.status,0,0,' ',0,0)

	elseif (MQH.name .eq. MQH_QUERY.name) then
	  if (MQD.QUERY.object .eq. MQH_BOXINFO.name) then
	    call INFOMSG (MSK_BOXES, info, sid)
	  elseif (MQD.QUERY.object .eq. MQH_COLORINFO.name) then
	    call INFOMSG (MSK_COLOR, info, sid)
	  elseif (MQD.QUERY.object .eq. MQH_MARKINFO.name) then
	    call INFOMSG (MSK_MARK, info, sid)
	  elseif (MQD.QUERY.object .eq. MQH_PANINFO.name) then
	    call INFOMSG (MSK_PANX, info+1, sid)
	  elseif (MQD.QUERY.object .eq. MQH_POINTINFO.name) then
	    call INFOMSG (MSK_POINTS, info, sid)
	  elseif (MQD.QUERY.object .eq. MQH_ZOOMINFO.name) then
	    call INFOMSG (MSK_ZOOM, info+1, sid)
	  elseif (MQD.QUERY.object .eq. '=SETVIEW') then
	    k = 1
	    i = Gx.view 
	    if (Gx.split) i = -i
	    MQD.dbuf(1) = Mx.stk(k).xmin
	    MQD.dbuf(2) = Mx.stk(k).xmax
	    MQD.dbuf(3) = Mx.stk(k).ymin
	    MQD.dbuf(4) = Mx.stk(k).ymax
	    MQD.dbuf(5) = Gx.plat
	    MQD.dbuf(6) = Gx.plon
	    MQD.dbuf(7) = Gx.pdis 
	    call M$SEND_MSGL ('=SETVIEW',0,sid,i,0,7,'D',MQD,0)
	  elseif (MQD.QUERY.object .eq. '=FILES') then
	    i = MQH.info
	    call XGEO_LIST (i,str,j,2)
	    MQD.sbuf(1:80) = str(1:max(1,j-1))
	    call M$SEND_MSGL ('=FILES',0,sid,i,0,1,'S[80]',MQD,0)
	  else
	    call M$WARNING ('Unknown QUERY object: '//MQD.QUERY.object)
	  endif

	elseif (MQH.name .eq. MQH_REFRESH.name) then
	  Gx.hold = .false.
	  call XGEO_REFRESH (0)

	elseif (MQH.name .eq. MQH_REREAD.name) then
	  k = MQH.info
	  call XGEO_RELOAD (P_LAYER,k)	  
	  call XGEO_REFRESH (0)

	elseif (MQH.name .eq. MQH_ZOOM.name) then
	  k = min(L_stack-1,Mx.level+1)
	  Mx.stk(k).xmin = MQD.ZOOM.xmin
	  Mx.stk(k).xmax = MQD.ZOOM.xmax
	  Mx.stk(k).ymin = MQD.ZOOM.ymin
	  Mx.stk(k).ymax = MQD.ZOOM.ymax
	  Mx.level = k
	  call XGEO_REFRESH (0)

	elseif (MQH.name .eq. MQH_UNZOOM.name) then
	  if (info.lt.0) info = max(0,Mx.level+info)
	  info = info+1 ! tranlate final zoom level to internal
	  do while (Mx.level.gt.info) 
	    k = Mx.level
	    if (Gx.bsmap(k).ip.ne.0) 
     &		call MX$MAP (Gx.bsmap(k),0,0,0,0,L_FreeMap)
	    Mx.level = k-1
	    if (Mx.level.eq.info) then
	      call XGEO_REFRESH (0)
	      call MSGDO (MSK_UNZOOM, Mx.level)
	    endif
	  enddo

	elseif (MQH.name .eq. MQH_PAN.name) then
	  k = Mx.level
	  Mx.stk(k).xmin = MQD.PAN.xmin
	  Mx.stk(k).xmax = MQD.PAN.xmax
	  Mx.stk(k).ymin = MQD.PAN.ymin
	  Mx.stk(k).ymax = MQD.PAN.ymax
	  call XGEO_REFRESH (1)

	elseif (MQH.name .eq. MQH_HARDCOPY.name) then
	  if (MX$HARDCOPY_OPEN(MQD.HARDCOPY.file, MQH.info)) then
	    call XGEO_REFRESH (-1)
	    call MX$HARDCOPY_CLOSE()
	  endif

	elseif (MQH.name .eq. MQH_NEWFILE.name) then
	  i = MQH.info
	  if (i.eq.0) i = 3
	  call XGEO_UNLOAD (i)
	  call XGEO_LOAD (i,MQD.NEWFILE.filename)
	  call XGEO_REFRESH (0)

	elseif (MQH.name .eq. MQH_OVERLAY.name) then
	  i = MQH.info
	  if (i.eq.0) i = 3
	  do n = 1,Gx.nlayer
	    ply = Gx.lloc(n)
	    if (ly.name .eq. MQD.OVERLAY.text) then
	      call XGEO_RELOAD (i,n)	  
	      call XGEO_REFRESH (0)
	      return
	    endif
	  enddo
	  call XGEO_LOAD (i,MQD.OVERLAY.text)
	  call XGEO_REFRESH (0)

	elseif (MQH.name .eq. MQH_SET.name .or. 
     &		MQH.name .eq. MQH_SETD.name) then
	  i = MQH.info
	  if (MQH.name .eq. MQH_SET.name) MQD.SETD.value = i
	  call M$MOVE (Gx.options,options,16) ! save current
	  call XGEO_SETD (MQD.SETD.object, i, MQD.SETD.value)
	  if (M$CHK(Gx.options,options,8).eq.0) Gx.bsmask = 0 ! changed BS
	  if (M$CHK(Gx.options,options,16).eq.0 .or.
     &	      MQH.name .eq. MQH_SETD.name) call XGEO_REFRESH (0)

	elseif (MQH.name .eq. MQH_TRACE.name) then
	  n = MQH.info
	  if (n.lt.0) n = Gx.nlayer +1+n
	  if (n.le.0 .or. n.gt.Gx.nlayer) call M$ERROR 
     &		('TRACE message to non-existent layer')
	  ply = Gx.lloc(n)
	  call MX$TRACE_MSG (MQD, ly.line, ly.thick, 
     &			ly.symbol, ly.radius, ly.color)
	  if (BTEST(ly.enable,PM_BS)) Gx.bsmask = 0
	  call XGEO_REFRESH (0)

	elseif (MQH.name .eq. MQH_T5REC.name .or.
     &		MQH.name .eq. MQH_T5TEXT.name) then
	  m = MQH.info
	  n = MQD.T5REC.layer
	  if (n.eq.0) then	! marker
	    pt5 = %loc( Gx.mrk )
	    i = 1
	  else
	    if (n.lt.0) n = Gx.nlayer +1+n
	    if (n.lt.0 .or. n.gt.Gx.nlayer) then
	call M$WARNING('T5REC/T5TEXT message to non-existent layer')
	      return
	    endif
	    ply = Gx.lloc(n)
	    if (ly.type .ne. 5000) then
	call M$WARNING('T5REC/T5TEXT message to non-type5000 layer')
	      return
	    endif
	    i = MQD.T5REC.index
	    if (i.lt.0) i = ly.nelem +1+i
	    if (i.lt.0 .or. i.gt.ly.nelem) then
	call M$WARNING('T5REC/T5TEXT message to non-existent index')
	      return
	    endif
	  endif
	  if (i.ne.0) then
	    i1 = i
	    i2 = i
	  else
	    i1 = 1
	    i2 = ly.nelem
	  endif
	  do i = i1,i2
	  if (n.ne.0) pt5 = ly.bufloc + (i-1) * L_t5recsz
	  if (BTEST(m,1)) t5.alt = MQD.T5REC.alt
	  if (BTEST(m,2)) t5.lat = MQD.T5REC.lat
	  if (BTEST(m,3)) t5.lon = MQD.T5REC.lon
	  if (BTEST(m,4)) t5.name = MQD.T5REC.name
	  if (BTEST(m,4)) t5.name_len = M$LENGTH(t5.name)
	  if (MQH.name .eq. MQH_T5REC.name) then
	    if (BTEST(m,5)) t5.azim = MQD.T5REC.azim
	    if (BTEST(m,6)) t5.elev = MQD.T5REC.elev
	    if (BTEST(m,7)) t5.roll = MQD.T5REC.roll
	    if (BTEST(m,8)) t5.r = MQD.T5REC.r
	    if (BTEST(m,9)) t5.a = MQD.T5REC.a
	    if (BTEST(m,10)) t5.b = MQD.T5REC.b
	    if (BTEST(m,11)) t5.color = MQD.T5REC.color
	    if (BTEST(m,12)) t5.symbol = MQD.T5REC.symbol
	    if (BTEST(m,13)) t5.size = MQD.T5REC.size
	    if (BTEST(m,14)) t5.cone = MQD.T5REC.cone
	    if (BTEST(m,15)) t5.invert = MQD.T5REC.invert
	    if (t5.r .eq.-1) call POSP_TO_AER (t5)
	  else
	    if (BTEST(m,5)) t5.text = MQD.T5TEXT.text
	    if (BTEST(m,5)) t5.text_len = M$LENGTH(t5.text)
	    if (BTEST(m,6)) t5.color = MQD.T5TEXT.color
	    if (BTEST(m,7)) t5.symbol = MQD.T5TEXT.symbol
	    if (BTEST(m,8)) t5.size = MQD.T5TEXT.size
	    t5.cone = P_TEXT
	    if (BTEST(m,10)) t5.invert = MQD.T5TEXT.invert
	  endif
	  enddo
	  if (BTEST(ly.enable,PM_BS)) Gx.bsmask = 0
	  call XGEO_REFRESH (0)

	elseif (MQH.name .eq. MQH_T5LAYER.name) then

	  n = MQH.info
	  if (n.lt.0) n = Gx.nlayer +1+n
	  if (MQD.T5LAYER.action .eq. 'CREATE') then
	    Gx.nlayer = Gx.nlayer+1
	    if (n.ne.0 .and. Gx.nlayer.ne.n) then
		call M$ERROR('Must create new layers in sequence')
		return
	    endif
	    n = Gx.nlayer
	    call M$MALLOC (L_layersz,Gx.lloc(n))
	    call XGEO_GET_DATA (%val(Gx.lloc(n)),' ',-2) ! to init layer
	  elseif (n.le.0 .or. n.gt.Gx.nlayer) then
	    call M$WARNING ('T5LAYER message to non-existent layer')
	    return
	  endif
	  ply = Gx.lloc(n)
	  i1 = MQD.T5LAYER.p1
	  i2 = MQD.T5LAYER.p2
	  if (MQD.T5LAYER.action .eq. 'DELETE') then
	    if (i1.eq.-1) i1 = ly.nelem
	    if (i2.eq.-1) i2 = ly.nelem
	    i1 = max(1,min(ly.nelem,i1))
	    i2 = max(i1,min(ly.nelem,i2))
	    call M$MOVE (%val(ly.bufloc+i2*ly.bpe),
     &		%val(ly.bufloc+(i1-1)*ly.bpe), (ly.nelem-i2)*ly.bpe)
	    ly.nelem = ly.nelem - (i2-i1+1)
	    if (BTEST(ly.enable,PM_BS)) Gx.bsmask = 0
	  elseif (MQD.T5LAYER.action .eq. 'COPY') then
	    if (i2.eq.-1) i2 = ly.nelem
	    if (i2.eq.0 .or. i2.gt.ly.nelem) then	! extending
	      i2 = ly.nelem+1
	      call M$REMALLOC_EXT (ly.bufsize,ly.bufloc,i2*ly.bpe)
	    endif
	    if (i1.gt.0 .and. i1.le.ly.nelem) then	! copy
	      call M$MOVE (%val(ly.bufloc+(i1-1)*ly.bpe),
     &		%val(ly.bufloc+(i2-1)*ly.bpe), ly.bpe)
	    else					! create
	      call T5INIT (%val(ly.bufloc+(i2-1)*ly.bpe), 
     &		P_TEXT, ly.color, ly.symbol, ly.radius)
	    endif
	    ly.nelem = max(i2,ly.nelem)
	  elseif (MQD.T5LAYER.action .eq. 'CREATE') then
	    ly.nelem = i1
	    ly.bufsize = ly.nelem*ly.bpe
	    call M$MALLOC (ly.bufsize, ly.bufloc)
	    do i = 1,ly.nelem
	      call T5INIT (%val(ly.bufloc+(i-1)*ly.bpe), 
     &		P_TEXT, ly.color, ly.symbol, ly.radius)
	    enddo
	  else
	    call M$ERROR ('Unsupported /T5LAYER/ message action: '//
     &				MQD.T5LAYER.action)
	  endif
	  call XGEO_REFRESH (0)

	elseif (MQH.name .eq. MQH_COLOR.name) then
	  k = MQH.info
	  if (k.ge.1 .and. k.le.Gx.ncolors) then
	    i = MQD.color.position
	    if (i.ge.0) call MX$COLORTAB (k, Mc.colormap(4), i)
	    if (i.lt.0) call MX$COLORTAB (k, MQD.COLOR, -2)
	  endif

	elseif (MQH.name .eq. 'SETVIEW') then
	  k = 1
	  Mx.level = k
	  Gx.view = iabs(MQH.info)
	  Gx.split = MQH.info .lt. 0
	  if (MQH.ndata .ge. 4) then
	    Mx.stk(k).xmin = MQD.dbuf(1)
	    Mx.stk(k).xmax = MQD.dbuf(2)
	    Mx.stk(k).ymin = MQD.dbuf(3)
	    Mx.stk(k).ymax = MQD.dbuf(4)
	  endif
	  if (MQH.ndata .ge. 7) then
	    Gx.plat = MQD.dbuf(5)
	    Gx.plon = MQD.dbuf(6)
	    if (MQD.dbuf(7).gt.1.01d0) Gx.pdis = MQD.dbuf(7)
	  endif
	  call XGEO_REFRESH (-1)

	elseif (MQH.name .eq. 'PIPESET') then
	  n = MQH.info
	  if (n.gt.0 .and. n.le.Gx.nlayer) then
	    ply = Gx.lloc(n)
	    lastmode = ly.pipemode    
	    if (ly.lhcb.ne.0) then
	      if (MQH.ndata.ge.1) ly.pipemode = MQD.dbuf(1)
	      if (MQH.ndata.ge.2) then
		if (MQD.dbuf(2).gt.0) then
		  ly.piperead = max(0,ly.piperead) + MQD.dbuf(2)
		else
		  ly.piperead = MQD.dbuf(2)
		endif
	      endif
	      if (MQH.ndata.ge.3) ly.pipelen = MQD.dbuf(3)
	      if (ly.pipemode.eq.PP_CLEAR) then
		ly.nelem = 0
		if (BTEST(ly.enable,PM_BS)) Gx.bsmask = 0
		call XGEO_REFRESH (0)
		ly.pipemode = lastmode
	      endif
	    endif
	  endif

	elseif (MQH.name .eq. 'MODTIME') then
	  Gx.time = MQD.dbuf(1)  		! new time
	  call XGEO_GET_SVS (%val(Gx.eloc(1)),Gx.time)
	  call M$DRSLT ('G_TIME',Gx.time)

	elseif (MQH.name .eq. 'MODFREQ') then
	  n1 = 1
	  n2 = Gx.nevent
	  if (MQH.info.gt.0 .and. MQH.info.le.Gx.nevent) then
	    n1 = MQH.info
	    n2 = MQH.info
	  elseif (MQH.info .ne. 0) then
	    call M$WARNING ('MODFREQ event # out of range')
	    n2 = 0
	  endif
	  do n = n1,n2
	    pev = Gx.eloc(n)
	    ev.flinkt = MQD.dbuf(1)
	    ev.flinka = MQD.dbuf(2)
	    ev.flinkb = MQD.dbuf(3)
	  enddo
	  call XGEO_REFRESH (0)

	elseif (MQH.name .eq. 'MARKER') then
	  Gx.mrk.y = MQD.dbuf(1)
	  Gx.mrk.x = MQD.dbuf(2)
	  Gx.mrk.z = MQD.dbuf(3)
	  if (MQH.ndata .ge. 4) then
	    Gx.mrk.name = MQD.abuf(4)
	    Gx.mrk.name_len = M$LENGTH(Gx.mrk.name)
	  endif
	  if (MQH.ndata .eq. 5) then
	    Gx.mrk.text = MQD.sbuf(32+1:32+48)		! avl dr 5112
	    Gx.mrk.text_len = M$LENGTH(Gx.mrk.text)
	  endif
	  if (MQH.info.ne.0) then
	    Gx.mrk.colorbg = MQH.info / 100
	    Gx.mrk.color = mod( MQH.info,100 )
	  endif
	  call XGEO_REFRESH (0)

	elseif (MQH.name .eq. 'AZEL') then
	  if (MQH.ndata .eq. 0) then
	    call XGEO_AZEL (info,Gx.mrk,MQD.dbuf(4))
	  else
	    call XGEO_AZEL (info,MQD.dbuf(1),MQD.dbuf(4))
	  endif
	  call M$SEND_MSGL ('=AZEL',0,sid,info,0,2,'D',MQD.dbuf(4),0)

	elseif (MQH.name .eq. 'GEOMODE') then
	  Gx.glmode = MQD.dbuf(1)
	  Gx.glevn1 = MQD.dbuf(2)
	  Gx.glevn2 = MQD.dbuf(3)
	  if (MQH.ndata.ge.4) Gx.glstd1 = MQD.dbuf(4)
	  if (MQH.ndata.ge.5) Gx.glstd2 = MQD.dbuf(5)

	elseif (MQH.name .eq. '?GEOLOC') then
	  pev = Gx.eloc(1)
	  Gx.time = MQD.dbuf(1)  ! new time
	  call XGEO_GET_SVS (ev,Gx.time)
	  ev.tdoa = MQD.dbuf(2)  ! tdoa		
	  ev.fdoa = MQD.dbuf(3)  ! fdoa		
	  Gx.glmode = 1	  
	  Gx.glevn1 = 1	  
	  Gx.glevn2 = 1	  
	  Gx.mrk.lat = MQD.dbuf(4)
	  Gx.mrk.lon = MQD.dbuf(5)
	  Gx.mrk.alt = MQD.dbuf(6)
	  k = XGEO_LOCATE ()
	  MQD.dbuf(4) = Gx.mrk.lat
	  MQD.dbuf(5) = Gx.mrk.lon
	  MQD.dbuf(6) = Gx.mrk.alt
	  call M$SEND_MSGL ('=GEOLOC',0,sid,k,0,6,'D',MQD,0)

	elseif (MQH.name .eq. '?LOCATE') then
	  k = XGEO_LOCATE ()
	  MQD.dbuf(1) = Gx.mrk.lat
	  MQD.dbuf(2) = Gx.mrk.lon
	  MQD.dbuf(3) = Gx.mrk.alt
	  call M$SEND_MSGL ('=LOCATE',0,sid,k,0,3,'D',MQD,0)

	elseif (MQH.name .eq. '?PREDICT') then
	  k = max(1,min(Gx.nevent,MQH.info))
	  pev = Gx.eloc(k)
	  call PREDICT_EVAL (0,ev,			!< event struct
     &	    MQD.dbuf(3),MQD.dbuf(1),MQD.dbuf(2),	!< alt,lat,lon
     &	    MQD.dbuf(1),MQD.dbuf(2),			!> tdoa,fdoa
     &	    MQD.dbuf(3),MQD.dbuf(4) )			!> sin(aoaa),sin(aoab) 
	  MQD.dbuf(3) = asin(MQD.dbuf(3))		! calc aoaa
	  MQD.dbuf(4) = asin(MQD.dbuf(4))		! calc aoab
	  call M$SEND_MSGL ('=PREDICT',0,sid,k,0,4,'D',MQD,0)

	elseif (MQH.name .eq. 'DUMPPRED') then
	  if (MQH.ndata.lt.1) MQD.dbuf(1) = Gx.tmin
	  if (MQH.ndata.lt.2) MQD.dbuf(2) = Gx.tmax
	  if (MQH.ndata.lt.3) MQD.dbuf(3) = 300.d0
	  if (MQH.ndata.lt.4) MQD.sbuf(24+1:24+80) = 'XGEO_PREDICTS'
	  i = max(1,min(Gx.nevent,MQH.info))
	  call XGEO_DUMP_PREDICTS (%val(Gx.eloc(i)), 
     &	    MQD.dbuf(1),MQD.dbuf(2),MQD.dbuf(3),MQD.sbuf(24+1:24+80))

	elseif (MQH.name .eq. 'EVENTMOD') then
	  k = iabs(MQH.info)
	  if (k.le.0 .or. k.gt.Gx.nevent) call M$ERROR
     &	  ('EVENTMOD message sent to non-existent event')
	  pev = Gx.eloc(k)
	  if (MQH.info .gt. 0) then
	    ev.enable = JIBSET (ev.enable,PM_GLOBAL)
	  else
	    ev.enable = JIBCLR (ev.enable,PM_GLOBAL)
	  endif
	  if (MQH.ndata.ge.1) ev.tdmode = MQD.dbuf(1)
	  if (MQH.ndata.ge.2) ev.fdmode = MQD.dbuf(2)
	  if (MQH.ndata.ge.3) ev.tdoa   = MQD.dbuf(3)
	  if (MQH.ndata.ge.4) ev.fdoa   = MQD.dbuf(4)
	  if (MQH.ndata.ge.5) ev.dtdoa  = MQD.dbuf(5)
	  if (MQH.ndata.ge.6) ev.dfdoa  = MQD.dbuf(6)
	  if (MQH.ndata.ge.7) ev.iballa = MQD.dbuf(7)
	  if (MQH.ndata.ge.8) ev.iballb = MQD.dbuf(8)
	  call XGEO_REFRESH (0)

	elseif (MQH.name .eq. 'EVENTREG') then
	  Gx.rlon1 = MQD.dbuf(1)
	  Gx.rlon2 = MQD.dbuf(2)
	  Gx.rlat1 = MQD.dbuf(3)
	  Gx.rlat2 = MQD.dbuf(4)
	  if (MQH.info .ne.-1) call XGEO_REFRESH (0)

	elseif (MQH.name .eq. MQH_ANNOTATE.name) then
	  call MX$FACILITIES (Gx.xf, 'A', 3)

	elseif (MQH.name(1:3).eq.'BOX' .or. MQH.name(1:5).eq.'POINT') then
	  if (MQH.name .eq. MQH_BOXDELETE.name .or.
     &	      MQH.name .eq. MQH_POINTDELETE.name) then
	    call MX$FACILITIES (Gx.xf, MQH.name(1:1), 4)
	  else
	    i = 2
	    if (MQH.info .le. -1) i = L_OneShot
	    if (MQH.name .eq. MQH_BOXADD.name) Gx.xf.iarg = Gx.xf.nbox+1
	    if (MQH.name .eq. MQH_POINTADD.name) Gx.xf.iarg = Gx.xf.npoint+1
	    if (MQH.name .eq. MQH_BOXADD.name .or.
     &		MQH.name .eq. MQH_POINTADD.name .or.
     &	      	MQH.name .eq. MQH_BOXALTER.name .or.
     &		MQH.name .eq. MQH_POINTALTER.name) 
     &		call MX$FACILITIES (Gx.xf, MQH.name(1:1), i)
	  endif

	elseif (MQH.name .eq. 'PAUSE') then
	  call M$PAUSE( float(MQH.info) )

	elseif (MQH.name .eq. 'DUMPEVF') then
	  call XGEO_EVENTS_SAVE ('XGEO_EVF',1)
	elseif (MQH.name .eq. 'LOADEVF') then
	  call XGEO_EVENTS_LOAD ('XGEO_EVF',-1)

	else
	  call M$WARNING('Unrecognized message: '//MQH.name)
	endif

	enddo

	if (M$RFIND('G_ACTION',action,dval) .eq. 'A') then
	  if (action .eq. 'DONE') return
	  if (action .eq. 'DUMPEVF') then
	    call XGEO_EVENTS_SAVE ('XGEO_EVF',1)
	  elseif (action .eq. 'LOADEVF') then
	    call XGEO_EVENTS_LOAD ('XGEO_EVF',-1)
	    call XGEO_REFRESH (0)
	  else
	    call M$WARNING('Unrecognized action result: '//action)
	  endif
	  Gx.status = Gx.status + 1
	  call M$ARSLT ('G_ACTION','DONE')
	endif

	return
	end


	subroutine MSGDO (bit, iarg)
	include 'xgeoplot.inc'
	integer*4 bit, iarg
	if (Gx.xf.msgid.ge.0 .and. BTEST(Gx.xf.msgmask,bit)) 
     &		call INFOMSG (bit, iarg, Gx.xf.msgid)
	return
	end

	subroutine INFOMSG (bit, iarg, rid)

	include 'xgeoplot.inc'
	include 'qmessdef.inc'
	include 'qmessdecl.inc'
	integer*4 iarg, rid, bit

	if (bit .eq. MSK_BOXES) then
	  Gx.xf.iarg = iarg
	  Gx.xf.mqdloc.p = %loc(MQD.BOXINFO)
	  call MX$FACILITIES (Gx.xf, 'B', 3) 
	  MQH_BOXINFO.boxid = Gx.xf.iarg
	  call M$SEND_MSG (MQH_BOXINFO, MQD.BOXINFO, rid, 0)
	elseif (bit .eq. MSK_POINTS) then
	  Gx.xf.iarg = iarg
	  Gx.xf.mqdloc.p = %loc(MQD.POINTINFO)
	  call MX$FACILITIES (Gx.xf, 'P', 3)
	  MQH_POINTINFO.pointid = Gx.xf.iarg
	  call M$SEND_MSG (MQH_POINTINFO, MQD.POINTINFO, rid, 0)
	elseif (bit .eq. MSK_COLOR) then
	  call MX$COLORTAB (iarg, MQD.COLORINFO, -3)
	  MQH_COLORINFO.colorid = iarg
	  call M$SEND_MSG (MQH_COLORINFO, MQD.COLORINFO, rid, 0)
	elseif (bit .eq. MSK_EXIT) then
	  call M$SEND_MSG (MQH_EXITINFO, MQD, rid, 0)
	elseif (bit .eq. MSK_MARK) then
	  MQD.MARKINFO.x = Gx.mrk.x
	  MQD.MARKINFO.y = Gx.mrk.y
	  MQD.MARKINFO.z = Gx.mrk.z
	  MQH_MARKINFO.info = Mx.state_mask
	  call M$SEND_MSG (MQH_MARKINFO, MQD.MARKINFO, rid, 0)
	elseif (bit .eq. MSK_ZOOM .or. bit .eq. MSK_UNZOOM) then
	  if (iarg.lt.1) i = max(1,Mx.level+iarg)
	  if (iarg.ge.1) i = min(iarg,Mx.level)
	  MQD.PLOTINFO.xmin = min (Mx.stk(i).xmin, Mx.stk(i).xmax)
	  MQD.PLOTINFO.xmax = max (Mx.stk(i).xmin, Mx.stk(i).xmax)
	  MQD.PLOTINFO.ymin = min (Mx.stk(i).ymin, Mx.stk(i).ymax)
	  MQD.PLOTINFO.ymax = max (Mx.stk(i).ymin, Mx.stk(i).ymax)
	  if (bit .eq. MSK_ZOOM) then
	    MQH_ZOOMINFO.level = i - 1
	    call M$SEND_MSG (MQH_ZOOMINFO, MQD.PLOTINFO, rid, 0)
	  elseif (bit .eq. MSK_UNZOOM) then
	    MQH_UNZOOMINFO.level = i - 1
	    call M$SEND_MSG (MQH_UNZOOMINFO, MQD.PLOTINFO, rid, 0)
	  endif
	else
	  call M$WARNING ('Unknown bit for INFOMSG')
	endif

	return
	end




	function XGEO_PAN_HANDLE (m,i) 
	include 'xgeoplot.inc'
	include 'xwidgets.inc'
	! m = 0:draw  1:check  2:drag  3:release
	! mode = 1:splitbox  2:vertical  3:horizontal  4:box  5/6:controldrag
	save
	real*8 d1,d2
	logical*4 act, inbox
	integer*4 junk,jx,jy,jb,panmode
	record /SCROLLBAR/ sbx,sby

c init returns
	i = 0
	k = Mx.level 
	act = .false.
	ix = Mx.xpos
	iy = Mx.ypos
	sbx.flag = -12
	sby.flag = -12

c handle calling mode
	if (m.eq.3) then
	  panmode = 0
	  goto 9
	elseif (m.eq.1) then
	  panmode = 0
	  sbx.action = 0
	  sby.action = 0
	  jb = Mx.button_press
	  if (jb.eq.0) goto 9
	  if (IAND(Mx.state_mask,L_ControlKeyMask).ne.0) jb = -jb
	elseif (m.eq.2) then
	  Mx.button_press = iabs(jb)
	endif
	if (.not. Gx.pan) goto 9
	if (k.eq.1) goto 7

c split box drag
	ix1 = ix - (Gx.zx+Gx.zw/2) 
	iy1 = iy - (Gx.zy+Gx.zh/2)
	if (m.eq.0) then
	  ! nothing to draw
	elseif (panmode.eq.1 .or. (panmode.eq.0 .and. 
     &		Gx.split .and. jb.eq.1 .and. 
     &		iabs(ix1).lt.10 .and. iabs(iy1).lt.10) )then
	  act = .true.
	  panmode = 1
	  if (iabs(ix1).gt.1 .or. iabs(iy1).gt.1) then
	    i = i+1
	    d1 = ix1 * Mx.stk(1).xscl
	    Mx.stk(k).xmin = Mx.stk(k).xmin + d1
	    Mx.stk(k).xmax = Mx.stk(k).xmax + d1
	    d1 = -iy1 * Mx.stk(1).yscl
	    Mx.stk(k).ymin = Mx.stk(k).ymin + d1
	    Mx.stk(k).ymax = Mx.stk(k).ymax + d1
	  endif
	endif

c vertical
	ix1 = Mx.stk(k).x2 + Mx.text_w
	ix2 = ix1 + Mx.text_w
	iy1 = Mx.stk(k).y1 + 3*Mx.text_w
	iy2 = Mx.stk(k).y2
	if (m.eq.0) then
	  junk = MX$SCROLLBAR (0,ix1,ix2,iy1,iy2,
     &	    Mx.stk(k).ymin,Mx.stk(k).ymax,Mx.stk(1).ymin,Mx.stk(1).ymax)
	elseif (panmode.eq.2 .or. (panmode.eq.0 .and.
     &		ix.gt.ix1-25*(m-1).and.iy.gt.iy1.and.iy.lt.iy2) )then
	  act = .true.
	  panmode = 2
	  d1 = Mx.stk(k).ymin
	  d2 = Mx.stk(k).ymax
	  junk = MX$SCROLLBAR (sby,ix1,ix2,iy1,iy2,
     &	    Mx.stk(k).ymin,Mx.stk(k).ymax,Mx.stk(1).ymin,Mx.stk(1).ymax)
	  if (Mx.stk(k).ymin.ne.d1 .or. Mx.stk(k).ymax.ne.d2) i = i+1
	endif

c horizontal
	ix1 = Mx.stk(k).x1
	ix2 = Mx.stk(k).x2 - 3*Mx.text_w
	if (Gx.axes) ix1 = min(ix2-40,ix1 + Mx.text_w*30)
	iy2 = Mx.stk(k).y1 - Mx.text_w
	iy1 = iy2 - Mx.text_w
	if (m.eq.0) then
	  junk = MX$SCROLLBAR (0,ix1,ix2,iy1,iy2,
     &	    Mx.stk(k).xmin,Mx.stk(k).xmax,Mx.stk(1).xmin,Mx.stk(1).xmax)
	elseif (panmode.eq.3 .or. (panmode.eq.0 .and.
     &		ix.gt.ix1.and.ix.lt.ix2.and.iy.lt.iy2+25*(m-1)) )then
	  act = .true.
	  panmode = 3
	  d1 = Mx.stk(k).xmin
	  d2 = Mx.stk(k).xmax
	  junk = MX$SCROLLBAR (sbx,ix1,ix2,iy1,iy2,
     &	    Mx.stk(k).xmin,Mx.stk(k).xmax,Mx.stk(1).xmin,Mx.stk(1).xmax)
	  if (Mx.stk(k).xmin.ne.d1 .or. Mx.stk(k).xmax.ne.d2) i = i+1
	endif

c box move
 7	ix1 = Mx.stk(k).x2 - 2*Mx.text_w
	iy1 = Mx.stk(k).y1 - 2*Mx.text_w
	if (k.eq.1 .and. .not. Gx.axes) then
	  ix1 = ix1 - 2*Mx.text_w
	  iy1 = iy1 + 2*Mx.text_w
	endif
	ix2 = ix1 + 4*Mx.text_w
	iy2 = iy1 + 4*Mx.text_w
	inbox = ix.ge.ix1.and.ix.le.ix2.and.iy.ge.iy1.and.iy.lt.iy2 
	if (k.eq.1 .and. (Gx.split .or. Gx.view.eq.1)) then
	  ! nobox
	elseif (m.le.0) then
	  call MX$TEXTLINE(ix1,iy1,ix2,iy1)
	  call MX$TEXTLINE(ix2,iy1,ix2,iy2)
	  call MX$TEXTLINE(ix2,iy2,ix1,iy2)
	  call MX$TEXTLINE(ix1,iy2,ix1,iy1)
	elseif ((panmode.eq.4 .or. panmode.eq.0) .and. inbox) then
	  act = .true.
	  panmode = 4
	endif

c control drag and box move handler
	if (panmode.ge.5 .or. (panmode.eq.4 .and. act)
     &	    .or. (panmode.eq.0 .and. jb.lt.0)) then
	  act = .true.
	  if (panmode.eq.0) then
	    panmode = 5
	    if (Gx.split .and. ix.le.Mx.stk(1).x2) k = 1
	    if (k.eq.1) panmode = 6
	    ! check for single click mode
	    fval = SECNDS(0.0)-.01
	    do while (SECNDS(fval).lt.0.5 .and.
     &		Mx.button_release .eq. 0 .and. 
     &		iabs(Mx.xpos-ix).le.5 .and. iabs(Mx.ypos-iy).le.5 )
	      if (.not. MX$IFEVENT() ) call M$PAUSE(Mc.pause)
	      if (Mx.button_release.ne.0) i = -2
	    enddo	
	    ix = Mx.xpos
	    iy = Mx.ypos
	    jx = ix
	    jy = iy
	  endif
 	  d1 = 0
	  d2 = 0
	  if (panmode.eq.6) k = 1
	  if (iabs(jb).eq.2 .or. (jb.eq.-1.and.Gx.view.eq.1)) then	! shift
	    if (i.eq.-2) then
	      call XGEO_SPECS ()
	      d1 = (ix-(Mx.stk(k).x1+Mx.stk(k).x2)/2) * Mx.stk(k).xscl
	      d2 = ((Mx.stk(k).y1+Mx.stk(k).y2)/2-iy) * Mx.stk(k).yscl
	    elseif (panmode.eq.4 .and. inbox) then
	      d1 = (2*ix-(ix1+ix2)) * Mx.stk(k).xscl
	      d2 = ((iy1+iy2)-2*iy) * Mx.stk(k).yscl
	      if (jb.lt.0) d1 = -d1
	      if (jb.lt.0) d2 = -d2
	    elseif (panmode.gt.4) then
	      d1 = (jx-ix) * Mx.stk(k).xscl
	      d2 = (iy-jy) * Mx.stk(k).yscl
	    endif
	    if (d1.ne.0.d0 .or. d2.ne.0.d0) i = i+1
	    Mx.stk(k).xmin = Mx.stk(k).xmin + d1
	    Mx.stk(k).xmax = Mx.stk(k).xmax + d1
	    Mx.stk(k).ymin = Mx.stk(k).ymin + d2
	    Mx.stk(k).ymax = Mx.stk(k).ymax + d2
	  elseif (jb.eq.-1 .and. Gx.view.gt.1) then	! rotate
	    if (i.eq.-2) then
	      call XGEO_SPECS ()
	      d1 = Gx.ret.lon - Gx.plon
	      d2 = Gx.ret.lat - Gx.plat
	    elseif (panmode.eq.4 .and. inbox) then
	      d1 = -20.0 * (2*ix-(ix1+ix2))/(8*Mx.text_w)
	      d2 = 20.0 * (2*iy-(iy1+iy2))/(8*Mx.text_w)
	    elseif (panmode.gt.4) then
	      ix1 = Mx.stk(k).x2-Mx.stk(k).x1
	      iy1 = Mx.stk(k).x2-Mx.stk(k).y1
	      d1 = 180.0 * (jx-ix) / min(ix1,iy1)
	      d2 = 180.0 * (iy-jy) / min(ix1,iy1)
	    endif
	    if (d1.ne.0.d0 .or. d2.ne.0.d0) i = i+1
	    Gx.plon = mod( Gx.plon+d1+540, 360.d0 ) - 180.d0
	    Gx.plat = mod( Gx.plat+d2+540, 360.d0 ) - 180.d0
	    k = 1
	  elseif (iabs(jb).eq.3 .or. jb.eq.1.and.inbox) then		! in/out
	    if (panmode.eq.4 .and. inbox) then
	      if (jb.lt.0) d2 = -0.1*(2*iy-(iy1+iy2))/(8*Mx.text_w)
	      if (jb.gt.0) d2 = 0.1 * dble(ix-ix1+iy-iy1)/(8*Mx.text_w)
	      if (jb.eq.1) d2 = -d2
	      i = i+1
	    elseif (panmode.gt.4) then
	      d2 = 2 * dble(jy-iy) / Mx.height
	      if (d2.ne.0.d0) i = i+1
	    endif
	    d1 = (Mx.stk(k).xmax - Mx.stk(k).xmin) * d2
	    Mx.stk(k).xmin = Mx.stk(k).xmin - d1
	    Mx.stk(k).xmax = Mx.stk(k).xmax + d1
	    d1 = (Mx.stk(k).ymax - Mx.stk(k).ymin) * d2 
	    Mx.stk(k).ymin = Mx.stk(k).ymin - d1
	    Mx.stk(k).ymax = Mx.stk(k).ymax + d1
	  endif
	  jx = ix
	  jy = iy
	  if (i.ne.0 .and. k.eq.1) Gx.bsmask = 0
	endif

 9	XGEO_PAN_HANDLE = act
	return
	end


	subroutine XGEO_END ()
	include 'xgeoplot.inc'
	call XGEO_UNLOAD (0)
	call XGEO_FACILITIES (-1)
	do k = 1,L_stack
	  if (Gx.bsmap(k).ip.ne.0) 
     &		call MX$MAP (Gx.bsmap(k),0,0,0,0,L_FreeMap)
	enddo
	i = M$DEL_RESULT ('G_ACTION',0)
	call MX$CLOSE
	return
	end


	subroutine XGEO_ZOOM (ix1,iy1,ix2,iy2)
	include 'xgeoplot.inc'
	real*8 ratio, d1, d2
	k = Mx.level
	if (Gx.split .and. k.gt.1 .and. ix1+ix2.lt.Mx.width) k=1
	if (k.eq.L_stack-1) then
	  call MX$MESSAGE ('No more stack levels',1.0)
	else
	  Mx.stk(k+1).xmin = Mx.stk(k).xmin + 
     &		(ix1-Mx.stk(k).x1)*Mx.stk(k).xscl
	  Mx.stk(k+1).xmax = Mx.stk(k).xmin + 
     &		(ix2-Mx.stk(k).x1)*Mx.stk(k).xscl
	  Mx.stk(k+1).ymin = Mx.stk(k).ymin + 
     &		(Mx.stk(k).y2-iy2)*Mx.stk(k).yscl
	  Mx.stk(k+1).ymax = Mx.stk(k).ymin + 
     &		(Mx.stk(k).y2-iy1)*Mx.stk(k).yscl
	  k = k+1
	  Mx.level = k
	  if (Gx.zfill .and. Gx.prop) then
	    call XGEO_REFRESH_VIEW (-9)
	    ratio = Mx.stk(k).xscl/Mx.stk(k).yscl  
	    if (ratio.gt.1) then
	      d1 = (Mx.stk(k).ymax+Mx.stk(k).ymin)/2
	      d2 = (Mx.stk(k).ymax-Mx.stk(k).ymin)/2
	      Mx.stk(k).ymin = d1 - d2*ratio
	      Mx.stk(k).ymax = d1 + d2*ratio
	    else
	      d1 = (Mx.stk(k).xmax+Mx.stk(k).xmin)/2
	      d2 = (Mx.stk(k).xmax-Mx.stk(k).xmin)/2
	      Mx.stk(k).xmin = d1 - d2/ratio
	      Mx.stk(k).xmax = d1 + d2/ratio
	    endif
	  endif
	  Gx.bsmap(k).mask = 0
	  call XGEO_REFRESH (1)
	  call MSGDO (MSK_ZOOM, Mx.level)
	endif
	return
	end



	subroutine XGEO_FACILITIES (n)
	include 'xgeoplot.inc'
	external XGEO_SPECS
c
c handle annotate, boxes and points facilities
	k = Mx.level
	if (n.eq.1 .or. n.eq.2) then	! draw
	  if (Gx.view.eq.1) then
	    if (Gx.xf.nbox.gt.0) call MX$FACILITIES (Gx.xf,'B',L_Draw)
	    if (Gx.xf.npoint.gt.0) call MX$FACILITIES (Gx.xf,'P',L_Draw)
	  endif
	  if (n.eq.1) then
	    if (Gx.xf.nannt.gt.0) call MX$FACILITIES (Gx.xf,'A',L_Draw)
	  endif
	elseif (n.eq.0) then	! init
#if _PGF77|_HPUX
	  Gx.xf.ROUTINE.p = get_address(XGEO_SPECS)
#else
	  Gx.xf.ROUTINE.p = %loc(XGEO_SPECS)
#endif
	  Gx.xf.xy = .true.
	  Gx.xf.mannt = 20
	  Gx.xf.mbox  = 20
	  Gx.xf.mpoint = 100
	  call MX$FACILITIES (Gx.xf,'BPA',L_InitFacilities)
	elseif (n.eq.-1) then	! exit
	  call MX$FACILITIES (Gx.xf,'BPA',L_CleanExit)
	elseif (n.eq.5) then	! 
	  call MX$FACILITIES (Gx.xf,'A',L_Interactive)
	elseif (Gx.view.ne.1) then
	  call MX$MESSAGE ('Points & Boxes in Lat/Lon only',-1.0)
	elseif (n.eq.3) then	! 
	  call MX$FACILITIES (Gx.xf,'B',L_Interactive)
	elseif (n.eq.4) then	! 
	  call MX$FACILITIES (Gx.xf,'P',L_Interactive)
	elseif (n.eq.-3) then	! 
	  call MX$FACILITIES (Gx.xf,'B',L_Oneshot)
	elseif (n.eq.-4) then	! 
	  call MX$FACILITIES (Gx.xf,'P',L_Oneshot)
	endif
	return
	end


	subroutine XGEO_SPECS ()
	include 'xgeoplot.inc'
	save
	character chara*116, charb*116, ca,cb, clab*9, rlab*9, date*20
	real*8 tdoa,fdoa,dx,dy,range,aoaa,aoab,azel(5)
	real*8 maoaa,maoab,saoaa,saoab
	integer*4 lc, lca,lct
	pointer (pev,ev)
c
c convert pixel coordinate to real-world alt/lat/lon
	call PIXTOGEO (Mx.xpos,Mx.ypos, Gx.ret.z, Gx.ret.y, Gx.ret.x)
c
c calculate alitiude
	if (Gx.lhadb.eq.0) then
	  Gx.ret.alt = 0
	else
	  call XGEO_ALTITUDE (%val(Gx.lhadb), Gx.ret)
	endif
c
c calculate range
	dx = Gx.ret.x-Gx.mrk.x
	dy = Gx.ret.y-Gx.mrk.y
	call GREATARC (Gx.ret,Gx.mrk,range)
	if (Gx.units .eq. P_ENGLISH) then
	  rlab = ' dR(sm): '
	  range = range / sm2meters
	elseif (Gx.units .eq. P_NAUTICAL) then
	  rlab = ' dR(nm): '
	  range = range / nm2meters
	else
	  rlab = ' dR(km): '
	  range = range / 1000
	endif
c
c update specs
	lc  = 54  ! length of noevent specs
	lca = 85  ! length of noAoA specs
	lct = lc  ! length for computing time insert
	tdoa = 0
	fdoa = 0
	aoaa = 0
	aoab = 0
	maoaa = 0
	maoab = 0
	saoaa = 0
	saoab = 0
	ca = CHAR( Gx.iballa + ICHAR('0') )
	cb = CHAR( Gx.iballb + ICHAR('0') )
	if (Gx.nevent .gt. 0) then
	  pev = Gx.eloc(1)
	  call PREDICT_EVAL_ECR (-1, ev, 
     &		ev.site(1),ev.site(2),ev.site(3),
     &		tdoa,fdoa,saoaa,saoab)
	  call PREDICT_EVAL (-1,ev,
     &		Gx.mrk.z,Gx.mrk.y,Gx.mrk.x,
     &		tdoa,fdoa,maoaa,maoab)
	  call PREDICT_EVAL (0,ev,
     &		Gx.ret.z,Gx.ret.y,Gx.ret.x,
     &		tdoa,fdoa,aoaa,aoab)
	  ca = CHAR( ev.iballa + ICHAR('0') )
	  cb = CHAR( ev.iballb + ICHAR('0') )
	  lc  = LEN(chara)
	  lct = lca
	  if (Gx.aoaspec) lct = lc
	endif
c
c display CEP or Altitude 
	if (Gx.cep) then
	  call CEP_EVAL (Gx.ret.alt,Gx.ret.lat,Gx.ret.lon,dval)
	  if (Gx.units.eq.P_ENGLISH) then
	    clab = ' CEP(sm):'
	    dval = dval / sm2meters
	  elseif (Gx.units.eq.P_NAUTICAL) then
	    clab = ' CEP(nm):'
	    dval = dval / nm2meters
	  else
	    clab = ' CEP(km):'
	    dval = dval / 1000
	  endif
	elseif (Gx.lhadb.ne.0) then
	  if (Gx.units.eq.P_ENGLISH.or.Gx.units.eq.P_NAUTICAL) then
	    clab = ' Alt(ft):'
	    dval = Gx.ret.alt / feet2meters
	  else
	    clab = ' Alt(m): '
	    dval = Gx.ret.alt
	  endif
	else
	  clab = ' Az(deg):'
	  call XGEO_AZEL (Gx.iballa,Gx.ret,azel)
	  dval = azel(1)
	endif
c
c display updated specs
	if (Gx.spec) then
	  ix = Mx.text_w
	  iy = Mx.height-Mx.text_h*1.2
	  write (chara,1001) Gx.ret.y,dy,clab,dval,
     &		ca,rad2deg*asin(aoaa),ca,cb,tdoa
	  write (charb,1002) Gx.ret.x,dx,rlab,range,
     &		cb,rad2deg*asin(aoab),ca,cb,fdoa
	  if (Gx.aoaspec) then
	    write (chara(lca+1:),1003) 
     &		ca,rad2deg*asin(saoaa),ca,rad2deg*asin(maoaa)
	    write (charb(lca+1:),1003) 
     &		cb,rad2deg*asin(saoab),cb,rad2deg*asin(maoab)
	  else
	    chara(lca+1:) = ' '	!so nothing leftover from last time
	    charb(lca+1:) = ' '
	  endif
 1001	  format('Lat:',f11.6,' dLat:',f11.6,a9,1pg12.6,
     &		' AoA',a1,':',0pf7.3,' Td',a1,a1,':',1pg13.6)
 1002	  format('Lon:',f11.6,' dLon:',f11.6,a9,1pg12.6,
     &		' AoA',a1,':',0pf7.3,' Fd',a1,a1,':',1pg13.6)
 1003	  format(2x,'SAoA',a1,':',0pf7.3,' MAoA',a1,':',0pf7.3)
	  if (Gx.units.eq.P_NAUTICAL) then
	    i = M$D2A_FORM (Gx.ret.y, chara(5:15), 'DMS')
	    i = M$D2A_FORM (Gx.ret.x, charb(5:15), 'DMS')
	    i = M$D2A_FORM (dy, chara(22:32), 'DMS')
	    i = M$D2A_FORM (dx, charb(22:32), 'DMS')
	  endif
	  call MX$TEXT(ix,iy,chara(1:lc))
	  iy = iy + Mx.text_h
	  call MX$TEXT(ix,iy,charb(1:lc))

	  ! if no room for the time specs - bail
	  if (ix+(lct+15)*Mx.text_w .gt. Mx.width) return

	  i = M$TIMES2STR (DNINT(Gx.time), 0.d0, ' ', date)
	  if (i.le.8) then
	    chara(1:9) = ' '
	    charb(1:9) = date(1:8)
	  elseif (i.lt.20) then
	    chara(1:9) = date(1:min(8,i-8))
	    charb(1:9) = date(i-7:i)
	  else
	    chara(1:9) = date(3:10)
	    charb(1:9) = date(13:20)
	  endif
	  ix = Mx.width - Mx.text_w*16
	  call MX$TEXT(ix,iy,' Time: '//charb(1:9))
	  iy = iy - Mx.text_h
	  call MX$TEXT(ix,iy,' Date: '//chara(1:9))
	endif
c
	return
	end


	subroutine PIXTOGEO (ix,iy, alt,lat,lon)
	include 'xgeoplot.inc'
	real*8 retx,rety,retz,rretx,rrety,rretz,p2
	real*8 rad,slat,clat,slon,clon,fact
	real*8 alt,lat,lon
c
c adjust/scale for segments and Y inversion
	k = Mx.level
	if (Gx.split .and. k.gt.1 .and. ix.lt.Mx.stk(1).x2) k=1
	retx = Mx.stk(k).xmin + (ix-Mx.stk(k).x1)*Mx.stk(k).xscl
	rety = Mx.stk(k).ymin + (Mx.stk(k).y2-iy)*Mx.stk(k).yscl
c
c viewpoint transformation
	if (Gx.view .eq.1) then
	  lon = retx
	  lat = min( 90.d0, max( -90.d0, rety ))
	else
	  p2 = Gx.pdis*Gx.pdis
	  slat = dsin(Gx.plat*ddeg2rad)
	  clat = dcos(Gx.plat*ddeg2rad)
	  rad = retx*retx + rety*rety
	  ! do inverse perspective stretch factor mapping
	  call QUADSOLVE (p2+rad, -2*p2, p2-1, fact)
	  if (fact.eq.0.d0) then	! off earth
	    rad = datan2 (rety,retx)
	    retx = cos(rad)
	    rrety = sin(rad)
	    rety =  rrety*clat
	    retz = -rrety*slat
	  else				! on earth
	    rrety = rety*fact
	    rretx = retx*fact
	    rretz = sqrt( 1.-rad*fact*fact)
	    retx = rretx
	    rety =  rrety*clat + rretz*slat
	    retz = -rrety*slat + rretz*clat
	  endif
	  lat = drad2deg*asin(rety)
	  retx = drad2deg*datan2(retx,retz)
	  lon = dmod( Gx.plon+retx+540.d0, 360.d0 ) - 180.d0
	endif
	alt = 0
	return
	end


	subroutine QUADSOLVE (a,b,c,x)
	real*8 a,b,c,x,d
	d = b*b - 4*a*c	
	if (d .le. 1.e-20 .or. abs(a) .lt. 1.e-20) then	! off earth
	  x = 0
	else			! on earth
	  x = ( -b - dsqrt(d) ) / (2*a)
	endif
	return
	end

	subroutine GREATARC (a,b,r)
	include 'constants.inc'
	real*8 a(3),b(3),r 
	real*8 ax,ay,az,ar,bx,by,bz,br,abr
	call GEO2ECR (a(1),a(2),a(3),ax,ay,az,ar)
	call GEO2ECR (b(1),b(2),b(3),bx,by,bz,br)
	r = 0
	if (ar*br .gt. 1.e-20)
     &	r = (ax*bx + ay*by + az*bz) / (ar*br)
	r = min(1.d0,max(-1.d0,r))
	r = acos(r) * earth_equatorial_radius
	return
	end

	subroutine XGEO_AZEL (iball,viewer,azel)
	include 'xgeoplot.inc'
	real*8 viewer(3), azel(2), quad(6), sv(6), r
	integer*4 iball
	azel(1) = 0
	azel(2) = 0
	if (iball.le.0 .or. iball.gt.Gx.nball) return
	quad(1) = viewer(1)	
	quad(2) = viewer(2)	
	quad(3) = viewer(3)	
	quad(4) = 0
	quad(5) = 0
	quad(6) = 0
	call XGEO_GET_SV (%val(Gx.lhsv(iball)), sv, Gx.time, Gx.site, 0)
	call M$ECR_TO_TOP (sv, 1, 1, quad)
	r = dsqrt( sv(1)*sv(1) + sv(2)*sv(2) + sv(3)*sv(3) )
	if (abs(sv(1)).gt.1e-20 .or. abs(sv(2)).gt.1e-20) then
	  azel(1) = datan2d( sv(2), sv(1) )			! azimuth
	else
	  azel(1) = 0
	endif
	if (r.gt.1e-20) then
	  azel(2) = dasind( min(1.d0,max(-1.d0,sv(3)/r)) )	! elevation
	else
	  azel(2) = 0
	endif
	return
	end


	subroutine XGEO_VEHICLE_CONE (flag)
	include 'xgeoplot.inc'
	record /LAYER/ lyv
	record /T5REC/ vt5(8)
c
	if (flag.eq.0) then
	  m = 0
	  do n = 1,Gx.nball
	   if (Gx.coneveh.eq.-1 .or. Gx.coneveh.eq.n) then
	    m = m+1
	    vt5(m).name = 'V'//CHAR(n+48)
	    vt5(m).name_len = 2
	    vt5(m).text_len = 0
	    call XGEO_GET_SV (%val(Gx.lhsv(n)),	vt5(m), Gx.time, Gx.site, 0)
	    call M$CAR_TO_GEO (vt5(m),1,1)
	    vt5(m).palt = Gx.mrk.alt
	    vt5(m).plat = Gx.mrk.lat
	    vt5(m).plon = Gx.mrk.lon
	    vt5(m).r = -1
	    call POSP_TO_AER (vt5(m))
	    vt5(m).a = Gx.coneang
	    vt5(m).b = Gx.coneang
	    if (Gx.svh) then
	      vt5(m).color = PC_HORIZON+(n-1)
	    else
	      vt5(m).color = PC_TRACE
	    endif
	    vt5(m).colorbg = 0
	    vt5(m).cone    = P_DEFCONE
	    vt5(m).invert  = 0
	    vt5(m).symbol  = L_xsymbol
	    vt5(m).size    = 2
	   endif
	  enddo
	  lyv.type = 5000
	  lyv.name = '*'
	  lyv.nelem = m
	  if (Gx.svh) then
	    lyv.color = PC_HORIZON+(n-1)
	  else
	    lyv.color = PC_TRACE
	  endif
	  lyv.tunits = 0
	  lyv.bufloc = %loc(vt5)
	  lyv.enable = P_ENABLE
	  call MX$DRAW_MODE (Gx.linew)
	  call XGEO_DRAW_LAYER (lyv)
	  call MX$DRAW_MODE (1)
	elseif (flag.eq.1) then

	endif
	return
	end



	subroutine GET_CONE (t5,cone,nlevels,npts)
	include 'xgeoplot.inc'
	record /T5REC/ t5
	record /POS/ cone(*)
	integer*4 npts,nlevels,conelev
	real*4 theta,dtheta,a2i,b2i
	real*4 ctheta,stheta,ctheta2,stheta2
	real*4 c,cr,cx,cy,cz,czc
	real*8 dfactor

	if (t5.cone.eq.P_NOCONE) then
	  return
	elseif (t5.cone.eq.P_DEFCONE) then
	  conelev = Gx.conelev
	else
	  conelev = t5.cone
	endif
	nlevels = iabs(conelev)

	npts = Gx.conediv+2
	dtheta = twopi/Gx.conediv
	if (abs(t5.a).gt.1.e-15) then
	  a2i = 1./(t5.a*t5.a)
	else
	  a2i = 1.e30
	endif
	if (abs(t5.b).gt.1.e-15) then
	  b2i = 1./(t5.b*t5.b)
	else
	  b2i = 1.e30
	endif
	cz = t5.r 
	czc = cz
	if (t5.invert) czc = 0

	cone(1).x = 0
	cone(1).y = 0
	cone(1).z = czc
	do i = 0,Gx.conediv
	  theta = i*dtheta
	  ctheta = cos(theta)
	  stheta = sin(theta)
	  ctheta2 = ctheta*ctheta
	  stheta2 = stheta*stheta
	  c = sqrt( 1./(a2i*ctheta2+b2i*stheta2) )
	  if (t5.r.eq.0) then		! ellipse
	    cr = c
	  else
	    cr = cz * tand(max(-89.0,min(89.0,c)))
	  endif
	  cx = cr * ctheta
	  cy = cr * stheta
	  cone(i+2).x = cx
	  cone(i+2).y = cy
	  cone(i+2).z = czc
	enddo
	do n = 2,nlevels
	  dfactor = DBLE(nlevels+1-n)/nlevels
	  czc = cz*dfactor
	  if (t5.invert) czc = cz - czc
	  nn = (n-1)*npts
	  do i = 1,npts
	    ii = i+nn
	    cone(ii).x = cone(i).x * dfactor
	    cone(ii).y = cone(i).y * dfactor
	    cone(ii).z = czc
	  enddo
	enddo
	ii = nlevels*npts+1
	cone(ii).x = 0
	cone(ii).y = 0
	cone(ii).z = 0
	if (t5.invert) cone(ii).z = cz
	call M$TOP_TO_ECR (cone,1,ii,t5)
	if (t5.r.ne.0) call CLIP_AT_EARTH (cone,ii,conelev)
	call M$CAR_TO_GEO (cone,1,ii)
	if (conelev.eq.P_GEOCONE) then
	  do i = 1,ii-1
	    if (cone(i).alt .gt. 0)  cone(i).alt = 0
	  enddo
	endif
	return
	end

	subroutine DVSMUL (a,b,c,n)
	real*8 a(*),b,c(*)
	do i = 1,n
	  c(i) = a(i)*b
	enddo
	return
	end

	subroutine CLIP_AT_EARTH (cone,nplot,conelev)

	include 'xgeoplot.inc'
	record /POS/ cone(*)
	integer*4 nplot,conelev
	real*4 xo,yo,zo,xd,yd,zd	! origin and delta vectors
	real*4 xe,ye,ze,xei,yei,zei	! eccentricity vector and inverse
	real*4 p,qA,qB,qC,qR,qS,qSR	! quadratic components
	real*4 poo,pod,pdd		! inner products

	xe = earth_equatorial_radius
	ye = earth_equatorial_radius
	ze = xe*sqrt(1.0-earth_eccentricity**2)
	xei = 1./xe
	yei = 1./ye
	zei = 1./ze
	xo = xei*cone(nplot).x
	yo = yei*cone(nplot).y
	zo = zei*cone(nplot).z
	do i = 1,nplot-1
	  xd = xei*cone(i).x - xo
	  yd = yei*cone(i).y - yo
	  zd = zei*cone(i).z - zo
	  pdd = xd*xd + yd*yd + zd*zd
	  pod = xo*xd + yo*yd + zo*zd
	  poo = xo*xo + yo*yo + zo*zo
	  qA = pdd
	  if (qA .lt. 1e-30) goto 10
	  qB = 2*pod
	  qC = poo - 1
	  qR = qB / (qA*2)	
	  qS = qC / qA
	  if (qR*qR .gt. qS) then
	    qSR = sqrt( qR*qR - qS )
	    p = -qR - qSR
	    if ((p.gt.0 .and. p.lt.1) .or. conelev.eq.P_GEOCONE) then
	      cone(i).x = (xo + p*xd) * xe
	      cone(i).y = (yo + p*yd) * ye
	      cone(i).z = (zo + p*zd) * ze
	    endif
	  elseif (conelev.eq.P_GEOCONE) then
	    qA = pod*pod - pdd
	    qB = 2*pod*(poo-1)
	    qC = poo*(poo-1)
	    qR = qB / (qA*2)
	    qS = qC / qA
	    qSR = sqrt( max(0.0,qR*qR - qS) )
	    if (qr.gt.0) then
	      p = -qR + qSR
	    else
	      p = -qR - qSR
	    endif
D	    write(*,*)  i, qr, qr*qr, qs, qsr, p
	    cone(i).x = (xo + p*xd) * xe
	    cone(i).y = (yo + p*yd) * ye
	    cone(i).z = (zo + p*zd) * ze
	  endif
  10	enddo
	return
	end

	subroutine XGEO_DRAW_LAYER (ly)
	include 'xgeoplot.inc'
c
	record /PIXSTRUCT/ pix(*),cpix(256)
	record /POS/ cone(1024)
	record /T5REC/ t5(1)
	pointer (pt5,t5),(ppix,pix)
	real*4	vreg,ttt
	integer*4 ic,jc,isy,isz,npts,nplot,nlevels
	logical*4 labels
c
	k = Mx.level
c	ttt = secnds(0.0)
c
c plot 5000 type data
	if (ly.type .eq. 5000) then
	  call M$MALLOC (ly.nelem*4, ppix)
	  pt5 = ly.bufloc 
	  call DRAW_VREG (vreg)
	  labels = Gx.labs .and. BTEST(ly.enable,PM_LABEL) .and.
     &		(ly.name.eq.'*' .or. 
     &	      ly.nelem .le. Gx.tlmax .or. vreg .le. Gx.tlreg)
	  nplot = ly.nelem
	  if (ly.tunits .eq. 1) then
	    call DRAWG (t5, 3,11, pix, nplot, 0, ly.color)
	  else
	    call DRAWG (t5, 3,11, pix, nplot, 0, 0)
	  endif
	  do j = 1,nplot
	    ic  = t5(j).color
	    jc  = t5(j).cone
	    isy = t5(j).symbol
	    isz = t5(j).size
	    if (ly.tunits.eq.1 .and. ly.tdelta*0.6 .lt.
     &	        abs( ly.tstart + (j-1)*ly.tdelta - Gx.time )) goto 10 
	    ii = pix(j).x
	    jj = pix(j).y
	    if (ii.ge.Mx.stk(k).x1 .and. ii.le.Mx.stk(k).x2 .and.
     &		jj.ge.Mx.stk(k).y1 .and. jj.le.Mx.stk(k).y2) then
 	      call MX$DRAW_SYMBOLS (ic,pix(j),1,isy,isz)
	      if (labels) then
		 if (t5(j).name_len.gt.0) call MX$COLORTEXT 
     &		   (ic, ii+5, jj-5, t5(j).name(1:t5(j).name_len) )
		i1 = 1
		kk = t5(j).text_len
		do while (jc.eq.P_TEXT .and. i1.le.kk)
		  i2 = INDEX(t5(j).text(i1:kk),'@') + i1-1
		  if (i2.eq.i1-1) i2 = kk + 1
		  call MX$COLORTEXT (ic, ii+5+(t5(j).name_len+1)
     &			*Mx.text_w, jj-5, t5(j).text(i1:i2-1) )
		  i1 = i2+1
		  jj = jj + Mx.text_h
		enddo
	      endif
	    endif
	    if (jc.ne.P_NOCONE .and. jc.ne.P_TEXT 
     &				.and. Gx.conemod.ne.0) then
	      call GET_CONE (t5(j),cone,nlevels,npts)
	      if (Gx.conemod .eq. 1 .or. Gx.conemod.eq.3) then
	       do n = 0,nlevels-1
		call DRAWG (cone(n*npts+2), 3,3, cpix, npts-1, 0, ic)
	       enddo
	      endif
cbug call GET_CONE again because something in above DRAWG wipes out CONE data
	      if (Gx.conemod.eq.3) call GET_CONE (t5(j),cone,nlevels,npts)
	      if (Gx.conemod .eq. 2 .or. Gx.conemod.eq.3) then
		call DRAWG (cone, 3,npts*3, cpix, nlevels+1, 0, ic)
	      endif
	    endif
 10	  enddo
c
c plot local 1999 type segments for data contours
	elseif (ly.type.eq.-1999) then
	  call M$MALLOC (L_pix*4, ppix)
	  call DRAWF (%val(ly.bufloc), 2,2, pix, 
     &		ly.nelem, ly.type, ly.color)
c
c plot line segment / polygon data 
	elseif (ly.type.eq.1999 .or. ly.type.eq.1998) then
	  call M$MALLOC (L_pix*4, ppix)
	  call DRAWF (%val(ly.bufloc), ly.dim,ly.spa, pix, 
     &		ly.nelem, ly.type, ly.color)
c
c plot 2000 type segments for terrain features
	elseif (ly.type.eq.2000) then
	  call M$MALLOC (ly.nsub*4, ppix)
	  call DRAWF2 (%val(ly.bufloc), pix, ly.nsub, ly.nelem, ly.type, 
     &		ly.color, ly.xstart, ly.xdelta, ly.ystart, ly.ydelta)

	else
	  call M$WARNING ('Unallowable layer type')
	endif
D	write(*,*) 'Time = ',secnds(ttt)

	call MX$XFLUSH()
	call M$FREE (ppix)
	return
	end



	subroutine DRAW_INIT ()
	include 'xgeoplot.inc'
c
	integer*4 L_vis
	parameter (L_vis=256)
	logical*4 kpix,kpixl
	record /PIXSTRUCT/ pix(*),ppix(4)
	byte   bbuf(*),btmp
	real*4 fbuf(*),ftmp
	real*8 dbuf(*),dtmp
	integer*4 nelem,istat,istatb,type,color
c
	SAVE
c
	integer*4 npts,nplot,ic,lat,lon,mode,itype,jdim,jspa
	real*4 slat,clat,slon,clon,mrkx,mrky,pdis,pang
	real*4 alon,alat,dlat,rlon,rlat,yvis,rvis
	real*4 xscale,yscale,Vtop,Vleft,vreg
	real*4 xmin,xmax,ymin,ymax,xoff,yoff
	real*4 x,y,z,xx,yy,zz,ff,f1,r,vr
	real*4 pplon,pplat,cp,sp,ct,ear,eari
	real*4 xstart,xdelta,ystart,ydelta
	integer*4 nvis,svis(L_vis),ivis(L_vis)
	integer*4 x1,x2,y1,y2,idmax,nx,ny
	integer*4 iix1,iix2,iiy1,iiy2,kx,ky,kxy
c
        real*4 tsin(-361:361),dtsin(-361:361)
        real*4 tcos(-361:361),dtcos(-361:361)
c
	mrkx = Gx.plon
	mrky = Gx.plat
	pdis = Gx.pdis
c
	k = Mx.level
	xmin   = Mx.stk(k).xmin
	xmax   = Mx.stk(k).xmax
	ymin   = Mx.stk(k).ymin
	ymax   = Mx.stk(k).ymax
	xscale = 1./Mx.stk(k).xscl
	yscale = -1./Mx.stk(k).yscl
	xoff   = 0
	yoff   = 0
	Vleft  = Mx.stk(k).xmin - Mx.stk(k).x1*Mx.stk(k).xscl
	Vtop   = Mx.stk(k).ymax + Mx.stk(k).y1*Mx.stk(k).yscl
	ix1 = Mx.stk(k).x1
	ix2 = Mx.stk(k).x2
	iy1 = Mx.stk(k).y1
	iy2 = Mx.stk(k).y2
	idmax = max( ix2-ix1, iy2-iy1 )/2
c
c perspective setup
	if (Gx.view .eq. 1) then
	  yvis = 0.
	  yy   = 1.
	  vreg = min(xmax-xmin,ymax-ymin)
	else
	  slat = sin(-mrky*deg2rad)
	  clat = cos(-mrky*deg2rad)
	  rlon = mrkx
	  pang = asin(1./pdis)
	  yvis = sin(pang)
	  rvis = (1./cos(pang))**2
	  vreg = 50*min(xmax-xmin,ymax-ymin)
	  ear  = earth_equatorial_radius
	  eari = 1./ear
	endif
	ff = 1
	return
c
c get vreg value
	entry DRAW_VREG (vr)
	vr = vreg
	return
c
c   create sin/cos tables
	entry SINCOS_INIT
        do i = -361,361
          tsin(i) = sin(float(i)*deg2rad)
          tcos(i) = cos(float(i)*deg2rad)
        enddo
        do i = -361,360
          dtsin(i) = tsin(i+1)-tsin(i)
          dtcos(i) = tcos(i+1)-tcos(i)
        enddo
	return

	entry DRAWF (fbuf,jdim,jspa,pix,nelem,type,color)
c
c plot segments
	j = 0
	mode = 0
	istatb = 0
	itype = iabs(type)
	if (type.eq.0) then
	  npts = nelem
	  mode = 1
	  j = -jspa*3
	endif
	do while (j.lt.nelem*jspa)
c
c prep segment mode 0=header 1=data
	if (mode.eq.0) then
	  npts  = fbuf(j+1)
	  nplot = 2
	  jj = j+jspa
	  if (fbuf(j+2).ne.-999 .and. istatb.eq.4) mode = 1
	  if (type.lt.0 .and. Gx.view.gt.1) mode = 1
	endif
	if (fbuf(j+2).eq.-998) then	! offset
	  xoff = fbuf(j+jspa+1)
	  yoff = fbuf(j+jspa+2)
	  if (Gx.view.eq.1) then
	    Vleft = Mx.stk(k).xmin-xoff - Mx.stk(k).x1*Mx.stk(k).xscl
	    Vtop  = Mx.stk(k).ymax-yoff + Mx.stk(k).y1*Mx.stk(k).yscl
	  endif
	  mode = 1
	  goto 90
	endif
	if (mode.eq.1) then
	  nplot = min(L_pix,npts)
	  jj = j+jspa*3
	  if (color.gt.0) then
	    ic = color
	  elseif (color.eq.PC_MAP) then
	    ic = fbuf(j+2) + PC_MAPOFF
	  elseif (color.eq.PC_USER) then
	    ic = fbuf(j+2) + PC_USEROFF
	  elseif (color.eq.PC_RAW) then
	    ic = fbuf(j+2)
	  else
	    ic = color
	  endif
	endif
c
c plot segment in Lat/Lon mode
	if (Gx.view .eq.1) then
	  if (mode.eq.0) then
	    istat = 0
	    if (fbuf(jj+1)+xoff .gt. xmax) goto 90
	    if (fbuf(jj+2)+yoff .gt. ymax) goto 90
	    if (fbuf(jj+jspa+1)+xoff .lt. xmin) goto 90
	    if (fbuf(jj+jspa+2)+yoff .lt. ymin) goto 90
	    istat = 1
	    if (fbuf(j+2).eq.-999 .and. 
     &		fbuf(jj+1)+xoff .gt. xmin .and.
     &		fbuf(jj+2)+yoff .gt. ymin .and.
     &		fbuf(jj+jspa+1)+xoff .lt. xmax .and.
     &		fbuf(jj+jspa+2)+yoff .lt. ymax ) istat = 4
	    goto 90
	  endif
	  do i = 1,nplot
	    ix = (fbuf(jj+1)-Vleft)*xscale
	    iy = (fbuf(jj+2)-Vtop)*yscale
	    pix(i).x = max(-16384,min(16384,ix))
	    pix(i).y = max(-16384,min(16384,iy))
	    jj = jj+jspa
	  enddo
	  if (color.eq.0) then
	  elseif (itype.eq.1999) then
	    call MX$DRAW_LINES (ic,pix,nplot)
	  elseif (itype.eq.1998) then
	    call MX$FILL_POLY (ic,pix,nplot,1)
	  endif
c
c plot segment in perspective mode
	else
	  if (mode.eq.0) then
	    x1 = ix2+100
	    y1 = iy2+100
	    x2 = ix1-100
	    y2 = iy1-100
	    istat = 0
	    if (nplot.eq.2) nplot = 4
	  endif
	  nvis = 0
	  kpixl = .false.
	  do i = 1,nplot
	    alon = fbuf(jj+1) - rlon + xoff
	    alat = fbuf(jj+2) + yoff
	    lat = alat
	    lon = alon
	    pplat = alat - lat
	    pplon = alon - lon
	    cp = tcos(lon) + pplon*dtcos(lon)
	    sp = tsin(lon) + pplon*dtsin(lon)
	    ct = tcos(lat) + pplat*dtcos(lat)
	    z  = tsin(lat) + pplat*dtsin(lat)
	    y = ct*cp
	    if (jdim.eq.2) then
	     yy = y*clat - z*slat
	     kpix = yy .gt. yvis
	     if (mode.eq.0 .and. .not.kpix) goto 10
	     xx = ct*sp
	     zz = y*slat + z*clat
	     ff = pdis/(pdis-yy)	! 2-D pdis-yy always > 0
	     xx = xx*ff
	     zz = zz*ff
	    else
	     r = (fbuf(jj+3)+ear)*eari
	     yy = r*(y*clat - z*slat)
	     xx = ct*sp
	     zz = y*slat + z*clat
	     kpix = yy.lt.pdis
	     if (kpix) ff = r*pdis/(pdis-yy)
	     xx = xx*ff
	     zz = zz*ff
	     kpix = kpix .and. (yy.gt.yvis .or. xx*xx+zz*zz.ge.rvis)
	     if (mode.eq.0 .and. .not.kpix) goto 10
	    endif
            ix = (xx-Vleft)*xscale
            iy = (zz-Vtop)*yscale
 10	    if (mode.eq.1) then	! data pass
              pix(i).x = max(-16384,min(16384,ix))
              pix(i).y = max(-16384,min(16384,iy))
	      if (kpix.neqv.kpixl) then
	       if (kpix) then
		nvis = nvis+1
		svis(nvis) = i
	       else
		ivis(nvis) = i - svis(nvis) 
	       endif
	      endif
 	      kpixl = kpix
	    else		! check pass
	      if (kpix) then
	        if (ix.lt.x1) x1 = ix
	        if (ix.gt.x2) x2 = ix
	        if (iy.lt.y1) y1 = iy
	        if (iy.gt.y2) y2 = iy
	        istat = istat + 1
	      endif
	      if (i.eq.2) then	! try opposite corners
	        jj = jj-jspa
		ftmp = fbuf(jj+2)
		fbuf(jj+2) = fbuf(jj+jspa+2)
		fbuf(jj+jspa+2) = ftmp
	        jj = jj-jspa
	      elseif (i.eq.4) then ! finish check
	        if (x1.gt.ix2 .or. y1.gt.iy2 .or.
     &		  x2.lt.ix1 .or. y2.lt.iy1) istat = 0
		if (fbuf(j+2).eq.-999 .and. istat.eq.4 .and.
     &	           (x1.lt.ix1 .or. y1.lt.iy1 .or.
     &		    x2.gt.ix2 .or. y2.gt.iy2) ) istat = 1
	        jj = jj-jspa
	        ftmp = fbuf(jj+2)
	        fbuf(jj+2) = fbuf(jj+jspa+2)
	        fbuf(jj+jspa+2) = ftmp
		goto 90
	      endif
	    endif
	    jj = jj+jspa
	  enddo
	  ! finish off last visual segment
	  if (kpix) ivis(nvis) = nplot+1 - svis(nvis) 
	  if (color.eq.0) then
	  elseif (itype.eq.1999) then
	    do i = 1,nvis
	      call MX$DRAW_LINES (ic,pix(svis(i)),ivis(i))
	    enddo
	  elseif (itype.eq.1998) then
	    if (nvis.eq.1) call MX$FILL_POLY (ic,pix,nplot,1)
	  endif
	endif
c
c finish the mode
 90	if (mode.eq.0) then
	  if (istat.eq.0) then
	    mode = 1
	  elseif (fbuf(j+2).eq.-999) then
	    istatb = istat
	    npts = 0
	    mode = 1
	  endif
	endif
	if (mode.eq.1) j = j + (3+npts)*jspa
	mode = 1-mode
	enddo
	! restore the offsets
	xoff   = 0
	yoff   = 0
	Vleft  = Mx.stk(k).xmin - Mx.stk(k).x1*Mx.stk(k).xscl
	Vtop   = Mx.stk(k).ymax + Mx.stk(k).y1*Mx.stk(k).yscl
	return


	entry DRAWG (dbuf,jdim,jspa,pix,nelem,type,color)

c
c plot segments
	j = 0
	mode = 0
	istatb = 0
	itype = iabs(type)
	if (type.eq.0) then
	  npts = nelem
	  mode = 1
	  j = -jspa*3
	endif
	do while (j.lt.nelem*jspa)
c
c prep segment mode 0=header 1=data
	if (mode.eq.0) then
	  npts  = dbuf(j+1)
	  nplot = 2
	  jj = j+jspa
	  if (dbuf(j+2).ne.-999 .and. istatb.eq.4) mode = 1
	  if (type.lt.0 .and. Gx.view.gt.1) mode = 1
	endif
	if (mode.eq.1) then
	  nplot = npts
	  jj = j+jspa*3
	  if (color.gt.0) then
	    ic = color
	  elseif (color.eq.PC_MAP) then
	    ic = dbuf(j+2) + PC_MAPOFF
	  elseif (color.eq.PC_USER) then
	    ic = dbuf(j+2) + PC_USEROFF
	  elseif (color.eq.PC_RAW) then
	    ic = dbuf(j+2)
	  else
	    ic = color
	  endif
	endif
c
c plot segment in Lat/Lon mode
	if (Gx.view .eq.1) then
	  if (mode.eq.0) then
	    istat = 0
	    if (dbuf(jj+3) .gt. xmax) goto 91
	    if (dbuf(jj+2) .gt. ymax) goto 91
	    if (dbuf(jj+jspa+3) .lt. xmin) goto 91
	    if (dbuf(jj+jspa+2) .lt. ymin) goto 91
	    istat = 1
	    if (dbuf(j+2).eq.-999 .and. 
     &		dbuf(jj+3) .gt. xmin .and.
     &		dbuf(jj+2) .gt. ymin .and.
     &		dbuf(jj+jspa+3) .lt. xmax .and.
     &		dbuf(jj+jspa+2) .lt. ymax ) istat = 4
	    goto 91
	  endif
	  nvis = 1
	  svis(nvis) = 1
	  do i = 1,nplot
	    ix = (dbuf(jj+3)-Vleft)*xscale
	    iy = (dbuf(jj+2)-Vtop)*yscale
	    if (i.gt.0 .and. (iabs(ix-pix(i-1).x).gt.idmax .or. 
     &		iabs(iy-pix(i-1).y).gt.idmax) .and. nvis.lt.L_vis) then
	      ivis(nvis) = i - svis(nvis) 
	      nvis = nvis+1
	      svis(nvis) = i
	    endif
	    pix(i).x = max(-16384,min(16384,ix))
	    pix(i).y = max(-16384,min(16384,iy))
	    jj = jj+jspa
	  enddo
	  ivis(nvis) = nplot+1 - svis(nvis) 
	  if (color.eq.0) then
	  elseif (itype.eq.0 .or. itype.eq.1999) then
	    do i = 1,nvis
	      call MX$DRAW_LINES (ic,pix(svis(i)),ivis(i))
	    enddo
	  elseif (itype.eq.1998) then
	    if (nvis.eq.1) call MX$FILL_POLY (ic,pix,nplot,1)
	  endif
c
c plot segment in perspective mode
	else
	  if (mode.eq.0) then
	    x1 = ix2+100
	    y1 = iy2+100
	    x2 = ix1-100
	    y2 = iy1-100
	    istat = 0
	    if (nplot.eq.2) nplot = 4
	  endif
	  nvis = 0
	  kpixl = .false.
	  do i = 1,nplot
	    r = (dbuf(jj+1)+ear)*eari
	    alat = dbuf(jj+2) 
	    alon = dbuf(jj+3) - rlon 
	    jj = jj+jspa
	    lat = alat
	    lon = alon
	    pplat = alat - lat
	    pplon = alon - lon
	    cp = tcos(lon) + pplon*dtcos(lon)
	    sp = tsin(lon) + pplon*dtsin(lon)
	    ct = tcos(lat) + pplat*dtcos(lat)
	    z  = tsin(lat) + pplat*dtsin(lat)
	    y = ct*cp
	    yy = r*(y*clat - z*slat)
	    kpix = yy.lt.pdis
	    xx = ct*sp
	    zz = y*slat + z*clat
	    if (kpix) ff = r*pdis/(pdis-yy)
	    xx = xx*ff
	    zz = zz*ff
	    kpix = kpix .and. (yy.gt.yvis .or. xx*xx+zz*zz.ge.rvis)
            ix = (xx-Vleft)*xscale
            iy = (zz-Vtop)*yscale
	    if (mode.eq.1) then	! data pass
	      if (.not.kpix) ix = -999
              pix(i).x = max(-16384,min(16384,ix))
              pix(i).y = max(-16384,min(16384,iy))
	      if (kpix.neqv.kpixl) then
	       if (kpix) then
		nvis = nvis+1
		svis(nvis) = i
	       else
		ivis(nvis) = i - svis(nvis) 
	       endif
	      endif
 	      kpixl = kpix
	    else		! check pass
	      if (ix.lt.x1) x1 = ix
	      if (ix.gt.x2) x2 = ix
	      if (iy.lt.y1) y1 = iy
	      if (iy.gt.y2) y2 = iy
	      istat = istat + 1
 	      if (i.eq.2) then	! try opposite corners
	        jj = jj-jspa-jspa
		dtmp = dbuf(jj+2)
		dbuf(jj+2) = dbuf(jj+jspa+2)
		dbuf(jj+jspa+2) = dtmp
	      elseif (i.eq.4) then ! finish check
	        if (x1.gt.ix2 .or. y1.gt.iy2 .or.
     &		  x2.lt.ix1 .or. y2.lt.iy1) istat = 0
		if (dbuf(j+2).eq.-999 .and. istat.eq.4 .and.
     &	           (x1.lt.ix1 .or. y1.lt.iy1 .or.
     &		    x2.gt.ix2 .or. y2.gt.iy2) ) istat = 1
	        jj = jj-jspa-jspa
	        dtmp = dbuf(jj+2)
	        dbuf(jj+2) = dbuf(jj+jspa+2)
	        dbuf(jj+jspa+2) = dtmp
		goto 91
	      endif
	    endif
	  enddo
	  ! finish off last visual segment
	  if (kpix) ivis(nvis) = nplot+1 - svis(nvis) 
	  if (color.eq.0) then
	  elseif (itype.eq.0 .or. itype.eq.1999) then
	    do i = 1,nvis
	      call MX$DRAW_LINES (ic,pix(svis(i)),ivis(i))
	    enddo
	  elseif (itype.eq.1998) then
	    if (nvis.eq.1) call MX$FILL_POLY (ic,pix,nplot,1)
	  endif
	endif
c
c finish the mode
 91	if (mode.eq.0) then
	  if (istat.eq.0) then
	    mode = 1
	  elseif (dbuf(j+2).eq.-999) then
	    istatb = istat
	    npts = 0
	    mode = 1
	  endif
	endif
	if (mode.eq.1) j = j + (3+npts)*jspa
	mode = 1-mode
	enddo
	return

c
c for drawing 2-D files (i.e. terrain maps)
	entry DRAWF2 (bbuf,pix,nx,ny,type,color,
     &		xstart,xdelta,ystart,ydelta)

	if (Gx.view .eq. 1) then
	  iix1 = max(1, INT((xmin-xstart)/xdelta)+1 )
	  iix2 = min(nx, INT((xmax-xstart)/xdelta)+2 )
	  if (ydelta.gt.0) then
	  iiy1 = max(1, INT((ymin-ystart)/ydelta)+1 )
	  iiy2 = min(ny, INT((ymax-ystart)/ydelta)+2 )
	  else
	  iiy1 = max(1, INT((ymax-ystart)/ydelta)+1 )
	  iiy2 = min(ny, INT((ymin-ystart)/ydelta)+2 )
	  endif
	  do j = iiy1, iiy2-1
	    x = xstart + (iix1-1)*xdelta 
	    ix = (x-Vleft)*xscale
	    ppix(2).x = ix
	    ppix(3).x = ix
	    y = ystart + (j-1)*ydelta
	    iy = (y-Vtop)*yscale
	    ppix(1).y = iy
	    ppix(2).y = iy
	    y = ystart + (j)*ydelta
	    iy = (y-Vtop)*yscale
	    ppix(3).y = iy
	    ppix(4).y = iy
	    ii = (j-1)*nx + iix1-1
	    do i = iix1, iix2-1
	      ppix(1) = ppix(2)
	      ppix(4) = ppix(3)
	      x = x + xdelta 
	      ix = (x-Vleft)*xscale
	      ppix(2).x = ix
	      ppix(3).x = ix
	      ii = ii+1
	      ic = PC_USEROFF + bbuf(ii)
	      call MX$FILL_POLY (ic,ppix,4,1)
	    enddo
	  enddo
	else
	 kxy = 8
	 do iiy1 = 1, ny-kxy, kxy
	 do iix1 = 1, nx-kxy, kxy
	  iiy2 = min(iiy1+kxy,ny)
	  iix2 = min(iix1+kxy,nx)
	  do mode = 0,1
	  if (mode.eq.0) then
	    kx = max(1,iix2-iix1)
	    ky = max(1,iiy2-iiy1)
	    x1 = ix2+100
	    y1 = iy2+100
	    x2 = ix1-100
	    y2 = iy1-100
	    istat = 0
	  else
	    if (istat.eq.0) goto 30
	    kx = 1
	    ky = 1
	  endif
	  do j = iiy1, iiy2, ky
	   ii = (j-1)*nx + iix1
	   do i = iix1, iix2, kx
	    alon = xstart + (i-1)*xdelta - rlon
	    alat = ystart + (j-1)*ydelta
	    lat = alat
	    lon = alon
	    pplat = alat - lat
	    pplon = alon - lon
	    cp = tcos(lon) + pplon*dtcos(lon)
	    sp = tsin(lon) + pplon*dtsin(lon)
	    ct = tcos(lat) + pplat*dtcos(lat)
	    z  = tsin(lat) + pplat*dtsin(lat)
	    y = ct*cp
	    yy = y*clat - z*slat
	    kpix = yy.gt.yvis
	    xx = ct*sp
	    zz = y*slat + z*clat
	    ff = pdis/(pdis-yy)		! 2-D (pdis-yy) always > 0 
	    xx = xx*ff
	    zz = zz*ff
            ix = (xx-Vleft)*xscale
            iy = (zz-Vtop)*yscale
	    if (mode.eq.0) then
	      if (ix.lt.x1) x1 = ix
	      if (ix.gt.x2) x2 = ix
	      if (iy.lt.y1) y1 = iy
	      if (iy.gt.y2) y2 = iy
	      if (kpix) istat = istat + 1
	    else
	      ppix(1) = ppix(2)
	      ppix(4) = ppix(3)
	      ppix(2) = pix(i)
	      ppix(3).x = ix
	      ppix(3).y = iy
	      if (kpix .and. i.gt.iix1 .and. j.gt.iiy1) then
	        ic = PC_USEROFF + bbuf(ii)
	        call MX$FILL_POLY (ic,ppix,4,1)
	      endif
	    endif
	    pix(i)  = ppix(3)
	    ii = ii+kx
	   enddo
	  enddo
	  if (mode.eq.0) then
	    if (x1.gt.ix2 .or. y1.gt.iy2 .or.
     &		x2.lt.ix1 .or. y2.lt.iy1) istat = 0
	  endif
  30	  enddo
	 enddo
	 enddo
	endif

	end


 	subroutine XGEO_GET_GRID (ly,array)
	include 'xgeoplot.inc'
	real*4 array(*)
c
	ly.nelem = L_grid
	ly.bufsize = L_grid*8
#if _SUN
	ly.bufloc = loc(array)
#else
	ly.bufloc = %loc(array)
#endif
	ly.type = -1999
	ly.color = PC_GRID
	ly.name = 'GRID'
c
	ii = 1
	do n = -75,75,15
	  array(ii) = 181			! # points
	  ii = ii+1
	  array(ii) = 0				! color ?
	  ii = ii+1
	  array(ii) = -180 + Gx.xoffs		! min lon
	  ii = ii+1
	  array(ii) = n				! min lat
	  ii = ii+1
	  array(ii) = 180 + Gx.xoffs		! max lon
	  ii = ii+1
	  array(ii) = n				! max lat
	  ii = ii+1
	  do nn = -180+Gx.xoffs,180+Gx.xoffs,2	! draw line at lat 'n'
	    array(ii) = nn
	    ii = ii+1
	    array(ii) = n
	    ii = ii+1
	  enddo
	enddo
	do n = -180+Gx.xoffs,165+Gx.xoffs,15
	  array(ii) = 89			! # points
	  ii = ii+1
	  array(ii) = 0				! color ?
	  ii = ii+1
	  array(ii) = n				! min lon
	  ii = ii+1
	  array(ii) = -88			! min lat
	  ii = ii+1
	  array(ii) = n				! max lon
	  ii = ii+1
	  array(ii) = 88			! max lat
	  ii = ii+1
	  do nn = -88,88,2			! draw line at lon 'n'
	    array(ii) = n
	    ii = ii+1
	    array(ii) = nn
	    ii = ii+1
	  enddo
	enddo

	ly.enable = P_ENABLE

	return
	end


	subroutine XGEO_GET_DATA (ly,name,mode)
	! mode: -1=init transform, 0=open, 1=reopen, 2=pipedata -2=init layer
	include 'xgeoplot.inc'
	include 'type5000.inc'
	character*(*) name
	record /HEADER/ hcb,hcbt,localhcb,hcbtx
	record /T5REC/ t5
	integer*4 mode, mask, color, ngot, nlose
	integer*4 q_name,q_aer,q_posp,q_elps,q_cone,q_flag,q_text
	real*4 zscale, zshift, zbuf(1)
	character*8 cmap
	byte bbuf(1)
	pointer (phcb,hcb),(pt5,t5),(pzbuf,zbuf),(pbbuf,bbuf)
	save hcbt
c
c branch for piped data read
	phcb = ly.lhcb
	if (mode.eq.2) goto 70
	phcb = %loc( localhcb )
c
c init transform data
	if (mode.eq.-1) then	
	  ! this is the template for the T5REC - if fields are not present
	  ! in the input files, they are initialized in M$TRANSFORM - so cool
	  call M$INITIALIZE (hcbt)
	  hcbt.type = 5000
	  hcbt.class = 5
	  hcbt.frame_of_ref = 'ECR'
	  hcbtx = hcbt
	  call M$ADDCOMP (hcbt,'POS ','VD',P_GEO,P_METERS)
	  call M$ADDCOMP (hcbt,'AER ','VD',P_CAR,P_DEGREES)
	  call M$ADDCOMP (hcbt,'CONE','VD',P_SPH,P_METERS)
	  call M$ADDCOMP (hcbt,'NAME','1A',0,0)
	  call M$ADDCOMP (hcbt,'FLAG','8B',0,0)
	  hcbt.bpe = M$BPE(hcbt)
	  call M$ADDCOMP (hcbtx,'POS ','VD',P_GEO,P_METERS)
	  call M$ADDCOMP (hcbtx,'TEXT','6A',0,0)
	  call M$ADDCOMP (hcbtx,'NAME','1A',0,0)
	  call M$ADDCOMP (hcbtx,'FLAG','8B',0,0)
	  return
	endif
c
c open input file
	call M$INIT (hcb,name,'1999,1998,2000,5000','S*,C*,V*,NH',0)
	call M$OPEN (hcb,HCBF_INPUT+HCBF_NOABORT)
	call TRIM_NAME (hcb.file_name, ly.name) 

	ly.type = 0
	ly.nelem = 0
	ly.bufloc = 0
	ly.tunits = hcb.tunits
	ly.tstart = hcb.tstart + hcb.timecode
	ly.tdelta = hcb.tdelta
	ly.piperead = 0
	ly.lhcb = 0
	ly.radius = 2
	ly.symbol = L_Xsymbol
	ly.line = 1
	ly.thick = 0
	ly.pipemode = PP_OFF
	ly.bpe = hcb.bpe

	if (mode.eq.-2) then
	  ly.type = 5000
	  ly.bpe = L_t5recsz
	  ly.color = PC_TRACE
	  ly.enable = IBCLR(P_DISABLE,PM_BS)
	  ly.name = 'CREATE'
	  return
	endif

	ly.enable = M$GET_QUALIFIER(hcb,'EM')
	if (ly.enable .eq. -1) then
	  ly.enable = P_ENABLE
	  ! if time varying or type 5000 - keep out of BS save
	  if (ly.tunits.eq.1 .or. hcb.class.eq.5) 
     &		ly.enable = IBCLR(ly.enable,PM_BS)
	endif	  

	if (.not. hcb.open) then
	  ly.nelem = -1

	elseif (hcb.type .eq. 1998 .or. hcb.type .eq. 1999) then
	  ly.type = hcb.type
	  ly.spa  = hcb.spa
	  ly.dim  = hcb.spa
	  ly.color = PC_MAP
	  if (M$GET_KEYWORD(hcb,'COLORMAP',cmap).gt.0) then
	    if (cmap.eq.'USER') ly.color = PC_USER
	    if (cmap.eq.'RAW') ly.color = PC_RAW
	  endif
	  ly.nelem = hcb.size
	  ly.bpe   = hcb.spa*4

	elseif (hcb.class .eq. 2) then
	  ly.type = 2000
	  ly.spa  = hcb.spa
	  ly.dim  = hcb.spa
	  ly.color = PC_USER
	  ly.xstart = hcb.xstart
	  ly.xdelta = hcb.xdelta
	  ly.ystart = hcb.ystart
	  ly.ydelta = hcb.ydelta
	  ly.nelem  = hcb.size
	  ly.nsub   = hcb.subsize
	  if (hcb.spa .ne. 1) call M$ERROR
     &	  ('Type 2000 data must be scalar colormap data scaled 1-16')
	  ly.bpe = hcb.subsize

	elseif (hcb.class .eq. 5) then
	  ly.type = 5000
	  ly.spa  = hcbt.components
	  ly.dim  = 3
	  if (mode.eq.0) ly.color = PC_TRACE
	  if (M$GET_KEYWORD(hcb,'COLORNAME',cmap) .gt. 0) then
	    if (cmap.eq.'GRID') ly.color = -1
	  endif
	  if (M$GET_KEYWORD(hcb,'COLORMAP',cmap).gt.0) then
	    if (cmap.eq.'USER') ly.color = PC_USER
	    if (cmap.eq.'RAW') ly.color = PC_RAW
	    if (cmap.eq.'MAP') ly.color = PC_MAP
	  endif
	  if (M$GET_KEYDATA(hcb,'COLOR',color,4,0,'L') .gt. 0) then
	    ly.color = color
	  endif
	  if (M$GET_KEYDATA(hcb,'THICK',color,4,0,'L') .gt. 0) then
	    ly.thick = color
	  endif
	  ly.nelem = hcb.size
	  ly.bpe  = L_t5recsz
	endif
	ly.bufsize = max(0, ly.nelem*ly.bpe)
	call M$MALLOC (ly.bufsize,ly.bufloc)
c
c check for pipe input and prep modes
	if (hcb.pipe .ne. 0) then
	  if (ly.nelem .lt. 0) call M$ERROR('Opening pipe: '//hcb.file_name)
	  call M$MALLOC (L_header_size, ly.lhcb)
	  call M$MOVE (hcb, %val(ly.lhcb), L_header_size)
	  ly.pipemode = PP_FIFO
	  ly.piperead = -1
	  ly.pipelen  = ly.nelem
	  ly.nelem = 0
	  return
	endif
c
c prepare for pipe input
  70	if (hcb.pipe .ne. 0) then
	  call M$HCBFUNC (hcb,HCBF_AVAIL,ngot)
	  if (ngot.le.0) return
	  if (ly.piperead.gt.0) ngot = min(ngot,ly.piperead)
	  if (ly.pipemode.eq.PP_WRAP .and. 
     &		ly.nelem.eq.ly.pipelen) ly.nelem = 0
	  if (ly.pipemode.eq.PP_FIFO .and. 
     &		ly.nelem+ngot.gt.ly.pipelen) then
	    nlose = max(0, min(ly.nelem, ly.nelem+ngot-ly.pipelen))
	    call M$MOVE (%val(ly.bufloc+nlose*ly.bpe),
     &		%val(ly.bufloc), (ly.nelem-nlose)*ly.bpe)
	    ly.nelem = ly.nelem - nlose
	  endif
	  if (ly.pipemode.eq.PP_CONSUME) then
	    call M$HCBFUNC (hcb,HCBF_SKIP,ngot)
	    return
	  endif
	  if (ly.pipelen.gt.0) ngot = min(ngot,ly.pipelen-ly.nelem)
	  hcb.xfer_len = ngot
	  call M$REMALLOC_EXT(ly.bufsize,ly.bufloc,(ly.nelem+ngot)*ly.bpe)
	else
	  hcb.xfer_len = ly.nelem
	  ly.nelem = 0
	endif
c
c read input data
	pbbuf = ly.bufloc + ly.nelem*ly.bpe
	ngot = 0
	if (.not. hcb.open) then
	  ngot = -1

	elseif (hcb.type .eq. 1998 .or. hcb.type .eq. 1999) then
	  call M$GRABX (hcb,bbuf,ngot)

	elseif (hcb.class .eq. 2) then
	  if (hcb.format .eq. 'SB') then
	    call M$GRABX (hcb,bbuf,ngot)
	  else
	    if (M$GET_KEYDATA(hcb,'ZSCALE',zscale,4,0,'F').ne.4) 
     &		zscale = 16.0 / (5000.0 + 5000.0)
	    if (M$GET_KEYDATA(hcb,'ZSHIFT',zshift,4,0,'F').ne.4) 
     &		zshift = 9.01
	    hcb.buf_type = 'F'
	    call M$MALLOC (hcb.xfer_len*max(hcb.bpe,hcb.ape*4),pzbuf)
	    call M$GRABX (hcb,zbuf,ngot)
	    do i = 1,ngot*hcb.ape
	      bbuf(i) = max(1.0, min(16.0, zbuf(i)*zscale+zshift ))
	    enddo
	    call M$FREE (pzbuf)
	  endif

	elseif (hcb.class .eq. 5) then
	  q_name = M$SUB_FIND (hcb,'NAME')
	  q_elps = M$SUB_FIND (hcb,'ELPS')
	  q_posp = M$SUB_FIND (hcb,'POSP')
	  q_aer  = M$SUB_FIND (hcb,'AER')
	  q_cone = M$SUB_FIND (hcb,'CONE')
	  q_flag = M$SUB_FIND (hcb,'FLAG')
	  q_text = M$SUB_FIND (hcb,'TEXT')
	  if (q_text.eq.0 .and. q_name.ne.0 .and. 
     &	      q_elps+q_posp+q_aer+q_cone.eq.0) then
	    if (M$NPT(hcb.comp(q_name).format) .gt. 1) then
	      hcb.comp(q_name).name = 'TEXT'
	      q_text = q_name
	      q_name = 0
	    endif
	  endif
	  if (q_text.gt.0) then
	    q_elps = 0
	    q_posp = 0
	    call M$GRABXT (hcb,bbuf,ngot,hcbtx)
	  else
	    if (q_posp.gt.0) then
	      hcbt.comp(2).name = 'POSP'
	      hcbt.comp(2).type = P_GEO
	      hcbt.comp(2).units = P_METERS
	    else
	      hcbt.comp(2).name = 'AER'
	      hcbt.comp(2).type = P_CAR
	      hcbt.comp(2).units = P_DEGREES
	    endif
	    if (q_elps.gt.0) then
	      hcbt.comp(3).name = 'ELPS'
	      hcbt.comp(3).type = 0
	      hcbt.comp(3).units = P_METERS
	    else
	      hcbt.comp(3).name = 'CONE'
	      hcbt.comp(3).type = P_SPH
	      hcbt.comp(3).units = P_METERS
	    endif
	    call M$GRABXT (hcb,bbuf,ngot,hcbt)
	  endif
	  do i = 1,ngot
	    pt5 = pbbuf + (i-1)*hcbt.bpe
	    if (q_flag.gt.0) then
	      ! t5.cone already set
	    else
	      t5.flags(1) = 0
	      t5.flags(2) = 0
	      if (q_cone.gt.0) t5.cone = P_DEFCONE 
	    endif
	    if (q_name.gt.0) then
	      if (t5.name_len.le.0) t5.name_len = M$LENGTH(t5.name)
	      if (t5.name_pos.le.0) t5.name_pos = 1
	    endif
	    if (q_text.gt.0) then
	      if (t5.text_len.le.0) t5.text_len = M$LENGTH(t5.text)
	      t5.cone = P_TEXT
	    endif
	    if (q_elps.gt.0) then
	      t5.cone = P_GEOCONE
	      t5.azim = 90-t5.tilt 	! in compass degrees
	      t5.elev = 0
	      t5.roll = 0
	      t5.b = t5.minor
	      t5.a = t5.major
	      t5.r = 0
	    endif
	    if (t5.color .gt. 0) then
	      if (ly.color .eq. PC_USER) then
		t5.color = t5.color + PC_USEROFF
	      elseif (ly.color .eq. PC_MAP) then
		t5.color = t5.color + PC_MAPOFF
	      endif
	    else
	      t5.color = ly.color
	    endif
	    if (t5.size.eq.0) t5.size = ly.radius
	    if (t5.symbol.eq.0) t5.symbol = ly.symbol
	    if (q_posp.gt.0) call POSP_TO_AER (t5)
D	    write(*,*) 'Pos = ',t5.alt,t5.lat,t5.lon
D	    write(*,*) 'AER = ',t5.azim,t5.elev,t5.roll
D	    write(*,*) 'CONE= ',t5.r,t5.a,t5.b
D	    write(*,*) 'LABL= ',t5.name,t5.name_len,t5.name_pos
D	    write(*,*) 'FLAG= ',t5.color,t5.symbol,t5.cone
	  enddo
	endif
	ly.nelem = ly.nelem + ngot
	if (ly.piperead.gt.0) ly.piperead = ly.piperead - ngot
c
c inform XGEO_PIPES of data update or close input file
	if (hcb.pipe.ne.0) then
	  Gx.status = Gx.status + 1
	else
	  call M$CLOSE (hcb)
	endif
	return

	end	  


	subroutine POSP_TO_AER (t5)
	include 'xgeoplot.inc'
	record /T5REC/ t5,t5a,t5b
	real*8 p(3),r,x,y,z

	t5a.alt = t5.alt
	t5a.lat = t5.lat
	t5a.lon = t5.lon
	t5a.azim = 0
	t5a.elev = 0
	t5a.roll = 0

	t5b.alt = t5.palt
	t5b.lat = t5.plat
	t5b.lon = t5.plon
	t5b.azim = 0
	t5b.elev = 0
	t5b.roll = 0

	p(1) = 0
	p(2) = 0
	p(3) = 0
	call M$TOP_TO_TOP (p,1,1,t5b,t5a)
	x = p(1)
	y = p(2)
	z = p(3)
	r = sqrt( x*x + y*y + z*z)
	t5.azim = 0
	if (x.ne.0 .or. z.ne.0) then
	  t5.elev = atan2d( -x, z )
	else
	  t5.elev = 0
	endif
	if (r.ne.0) then
	  t5.roll = asind( min(1.d0,max(-1.d0, -y/r )))
	else
	  t5.roll = 0
	endif
	if (t5.r .eq. -1) t5.r = r

	return
	end

	subroutine XGEO_LAYER_MOD (k)
	include 'xgeoplot.inc'
	pointer (ply,ly)

	ply = Gx.lloc(k)
	m = MX$MENU(0,'LAYER;Enables,Delete,Rename,Color')
	if (m.eq.1) then
	  if (XGEO_ENABLED(ly.enable,-9)) continue
	elseif (m.eq.2 .and. k.eq.0) then
	  call MX$MESSAGE ('GRID layer is special - cannot delete',-1.0)
	elseif (m.eq.2) then
	  Gx.bsmask = 0
	  call M$FREE (ly.bufloc)
	  call M$FREE (Gx.lloc(k))
	  Gx.nlayer = Gx.nlayer - 1
	  do k = k,Gx.nlayer
	    Gx.lloc(k) = Gx.lloc(k+1)
	  enddo
	elseif (m.eq.3) then
	  call MX$PROMPT ('New Name: ',ly.name)
	elseif (m.eq.4) then
	  m = ly.color
	  call MX$LVAL ('Color: ',ly.color,-1,Gx.ncolors,1,0)
	  if (ly.color.ne.m) Gx.bsmask = 0
	endif
	end


	subroutine XGEO_EVENT_INIT (k)
	include 'xgeoplot.inc'
	pointer (pev,ev)

	pev = Gx.eloc(k)

	ev.name = 'ROOT'
	ev.gsmode = PG_CURRENT
	ev.svmode = PG_CURRENT
	ev.tdmode = PG_CONTOUR
	ev.fdmode = PG_CONTOUR
	ev.iballa = Gx.iballa
	ev.iballb = Gx.iballb
	ev.enable = P_ENABLE
	ev.flags  = -1
	ev.dtime  = 0
	ev.tdoa   = 0
	ev.fdoa   = 0
	ev.dtdoa  = 0
	ev.dfdoa  = 0
	if (M$RFIND('G_FLINKT',str,ev.flinkt).le.'A') ev.flinkt = 100.d6
	if (M$RFIND('G_FLINKA',str,ev.flinka).le.'A') ev.flinka = -100.d6
	if (M$RFIND('G_FLINKB',str,ev.flinkb).le.'A') ev.flinkb = -100.d6
	call M$MOVE (Gx.site, ev.site, 24)
	call XGEO_GET_SVS (ev, Gx.time)

	return
	end


	subroutine XGEO_EVENT_MOD (k)
	include 'xgeoplot.inc'
	pointer (pev,ev)

	pev = Gx.eloc(k)
	m = MX$MENU(0,'EVENT;Enables,Delete,Modify')
	if (m .eq. 1) then
	  if (XGEO_ENABLED(ev.enable,-9)) continue
	elseif (m.eq.2) then
	  if (Gx.nevent .gt. 1) then
	    call M$FREE (Gx.eloc(k))
	    call M$MOVE (Gx.eloc(k+1),Gx.eloc(k),(Gx.nevent-k)*4)
	    Gx.nevent = Gx.nevent-1
	    if (k.eq.1) then		! root event
	      pev = Gx.eloc(k)
	      Gx.iballa = ev.iballa
	      Gx.iballb = ev.iballb
	      call XGEO_GET_SVS (ev,Gx.time)
	    endif
	  else
	    call MX$MESSAGE ('Must leave one event - not deleting',-1.0)
	  endif
	  return
	elseif (m.eq.3) then
	  call XGEO_EVENT_DIALOG (k)
	endif
	return
	end


	subroutine XGEO_EVENT_DIALOG (k)
	include 'xgeoplot.inc'

	if (Gx.nevent.le.0) then
	  call MX$MESSAGE ('No events loaded',1.0)
	else
	  call XGEO_EVENTS_SAVE ('XGEO_EVF',k)
	  call M$DISPATCH ('XGEOEVENT XGEO_EVF',2,-1)
	endif

	return
	end



	subroutine XGEO_GET_TIME (prompt,time)
	real*8 time
	character str*32, prompt*(*)
	str = ' '
	i = M$SEC2TOD(time,str)
	call MX$PROMPT (prompt,str)
	i = M$LENGTH(str)
	j = M$TOD2SEC(str(1:i),time)
	return
	end


	subroutine XGEO_DRAW_EVENT (n)
	include 'xgeoplot.inc'
	real*8 tdoa(3),fdoa(3),maoaa,maoab,saoaa,saoab
	pointer (pev,ev)
	l = Gx.lmg
	m = Gx.nmg*Gx.nmg*4
	pev = Gx.eloc(n)
	if (Gx.sview) 
     &	  call PREDICT_EVAL_ECR (-1,ev,ev.site(1),ev.site(2),ev.site(3),
     &		tdoa,fdoa,saoaa,saoab)
	if (Gx.mview) 
     &	  call PREDICT_EVAL (-1,ev,Gx.mrk.z,Gx.mrk.y,Gx.mrk.x,
     &		tdoa,fdoa,maoaa,maoab)
	if ((.not.Gx.mview .or. (maoaa.ge.0 .and. maoab.ge.0)) .and. 
     &	    (.not.Gx.sview .or. (saoaa.ge.0 .and. saoab.ge.0))) 
     &	  call XGEO_DRAW_EVENT_ (n,ev,Gx.nmg,
     &		%val(l+m*0),%val(l+m*1),%val(l+m*2),
     &		%val(l+m*3),%val(l+m*4))
	return
	end


	subroutine XGEO_DRAW_EVENT_ (n,ev,mg,td,fd,aa,ab,pchk)
	include 'xgeoplot.inc'
	integer*4 lne,lnf,lnt,mg,repeat
	real*4 td(*),fd(*),aa(*),ab(*),pchk(*)
	real*4 tmin,tmax,fmin,fmax
	real*8 lat,lon,alt,tdoa,fdoa,aoaa,aoab,tdoam,fdoam
	real*8 xs,xd,ys,yd, x,y,z
	real*4 cs,cd,cxs,cxd,cys,cyd
	real*4 buf(L_ctmp/4+6)
	character stmp
c
c calculate the grid boundaries
	k = Mx.level
	if (Gx.view.eq.1) then	! lat/lon
	  xs = Mx.stk(k).xmin
	  xd = Mx.stk(k).xmax 
	  ys = Mx.stk(k).ymin 
	  yd = Mx.stk(k).ymax 
	else			! perspective
	 if (k.gt.1) then
	  call PIXTOGEO (Mx.stk(k).x1,Mx.stk(k).y1, z,yd,xd)
	  xs = xd
	  ys = yd	
	  call PIXTOGEO (Mx.stk(k).x2,Mx.stk(k).y1, z,y,x)
	  xs = min(x,xs)
	  xd = max(x,xd)
	  ys = min(y,ys)
	  yd = max(y,yd)
	  call PIXTOGEO (Mx.stk(k).x2,Mx.stk(k).y2, z,y,x)
	  xs = min(x,xs)
	  xd = max(x,xd)
	  ys = min(y,ys)
	  yd = max(y,yd)
	  call PIXTOGEO (Mx.stk(k).x1,Mx.stk(k).y2, z,y,x)
	  xs = min(x,xs)
	  xd = max(x,xd)
	  ys = min(y,ys)
	  yd = max(y,yd)
	 endif
	 if (k.eq.1 .or. xd-xs.gt.180) then	! got a pole
	  xs = -180
	  xd = +180
	  y = 0.5*(ys+yd)
	  if (k.eq.1 .or. y.gt.0) yd = +90
	  if (k.eq.1 .or. y.lt.0) ys = -90
	 endif
	endif
	if (Gx.rlon1 .gt. xs) xs = Gx.rlon1
	if (Gx.rlon2 .lt. xd) xd = Gx.rlon2
	if (Gx.rlat1 .gt. ys) ys = Gx.rlat1
	if (Gx.rlat2 .lt. yd) yd = Gx.rlat2
	ys = max(-89.d0,ys)	! watch the poles
	yd = min(+89.d0,yd)	! watch the poles
D	write(*,*)  'Region = ',xs,xd,ys,yd
	xd = (xd - xs) / (mg-1)
	yd = (yd - ys) / (mg-1)
c
c set up the grid offsets for precision plotting real*4 variables
	cxs = 0
	cys = 0
	cxd = xd
	cyd = yd
	buf(1) = 0	! # points in segment
	buf(2) = -998	! special code for specifying offset
	buf(3) = xs	! bottom-left of grid
	buf(4) = ys
	alt = 0		! set altitude of contours to zero
	lnt = 0
	lnf = 0
	ly.type = -1999
	ly.name = 'EVENT'
	ly.bufloc = %loc(buf)
c
c set line width
	call MX$DRAW_MODE (Gx.linew)
c
c branch for computing CEP contours
	if (n.le.0) then
	  fmin = +1e30
	  fmax = -1e30
	  ii = 0
	  do j = 1,mg
	    lat = ys + yd*(j-1)
	    do i = 1,mg
	      lon = xs + xd*(i-1)
	      ii = ii+1
	      call CEP_EVAL  (alt,lat,lon, dval)
	      fd(ii) = dval
	    enddo
	  enddo
	  cs = 0
	  if (n.eq.0) then
	    cd = Gx.cepstep
	    nn = 30
	  endif
	  call MX$CONTOUR (cs,cd,nn,fd,pchk,
     &      buf(7), lne, L_ctmp, mg,mg,cxs,cxd,cys,cyd)
	  if (lne.gt.0) then
	    ly.nelem = lne+3
	    ly.color = PC_CEP
	    ly.enable = P_ENABLE
	    call XGEO_DRAW_LAYER (ly)
	  endif
	  return
	endif
c
c copy event enable flags to layer enable
	ly.enable = ev.enable
c
c get new TDOA/FDOA ?
	if (BTEST(ev.flags,PM_AUTOTF))  call PREDICT_EVAL (0,ev, 
     &	  Gx.mrk.alt,Gx.mrk.lat,Gx.mrk.lon, ev.tdoa,ev.fdoa, aoaa,aoab)
c
c compute t/fdoa offset for data grids in contour modes
  1	repeat = 0
	lat = ys + yd*mg/2
	lon = xs + xd*mg/2
	if (ev.tdmode .eq. PG_CONTOUR .and.  ev.tdoa.eq.0) then
	  call PREDICT_EVAL (0,ev, alt,lat,lon, tdoam,dval, aoaa,aoab)
	  if (ev.dtdoa.ne.0) tdoam = DNINT(tdoam/ev.dtdoa)*ev.dtdoa
	else
	  tdoam = ev.tdoa
	endif
	if (ev.fdmode .eq. PG_CONTOUR .and.  ev.fdoa.eq.0) then
	  call PREDICT_EVAL (0,ev, alt,lat,lon, dval,fdoam, aoaa,aoab)
	  if (ev.dfdoa.ne.0) fdoam = DNINT(fdoam/ev.dfdoa)*ev.dfdoa
	else
	  fdoam = ev.fdoa
	endif
c
c fill in the data grids
	tmin = +1e30
	tmax = -1e30
	fmin = +1e30
	fmax = -1e30
	ii = 0
	do j = 1,mg
	  lat = ys + yd*(j-1)
	  do i = 1,mg
	    lon = xs + xd*(i-1)
	    ii = ii+1
	    call PREDICT_EVAL (0,ev, alt,lat,lon, tdoa,fdoa, aoaa,aoab)
	    td(ii) = tdoa-tdoam
	    fd(ii) = fdoa-fdoam
	    aa(ii) = aoaa
	    ab(ii) = aoab
	    if (td(ii).lt.tmin) then
	      tmin = td(ii)
	    elseif (td(ii).gt.tmax) then
	      tmax = td(ii)
	    endif
	    if (fd(ii).lt.fmin) then
	      fmin = fd(ii)
	    elseif (fd(ii).gt.fmax) then
	      fmax = fd(ii)
	    endif
	  enddo
	enddo
c
c catch the zero dtfdoa cases
	if (ev.tdmode .eq. PG_CONTOUR .and. ev.dtdoa.eq.0) then
	  ev.dtdoa = (tmax-tmin)/20
	  repeat = repeat + 1
	endif
	if (ev.fdmode .eq. PG_CONTOUR .and. ev.dfdoa.eq.0) then
	  ev.dfdoa = (fmax-fmin)/20
	  repeat = repeat + 1
	endif
	if (repeat .gt. 0) goto 1
c
c draw the TDOA line,region,contour as calculated
	if (ev.tdmode .eq. PG_CONTOUR) then
	  if (ev.dtdoa .gt. 0.d0) then	! use delta
	    cd = ev.dtdoa
	  elseif (ev.dtdoa .lt. 0.d0) then	! use # contours
	    cd = (tmin-tmax)/ev.dtdoa
	  endif
	  cs = INT(tmin/cd)*cd
	  nn = (tmax-cs)/cd + 1
	elseif (ev.tdmode .eq. PG_REGION) then
	  cs = -ev.dtdoa/2
	  cd = ev.dtdoa
	  nn = 2
	elseif (ev.tdmode .eq. PG_LINE) then
	  cs = 0
	  cd = ev.dtdoa
	  nn = 1
	else
	  nn = 0
	endif
	if (nn.gt.0) then
	  call MX$CONTOUR (cs,cd,nn,td,pchk,
     &      buf(7), lnt, L_ctmp, mg,mg,cxs,cxd,cys,cyd)
	  if (lnt.gt.0) then
	    ly.nelem = lnt+3
	    ly.color = PC_TDOA
	    call XGEO_DRAW_LAYER (ly)
	    Gx.evp(4).lon = buf(13) + xs
	    Gx.evp(4).lat = buf(14) + ys
	  endif
	endif
c
c draw the FDOA line,region,contour as calculated
	if (ev.fdmode .eq. PG_CONTOUR) then
	  if (ev.dfdoa .gt. 0.d0) then	! use delta
	    cd = ev.dfdoa
	  elseif (ev.dfdoa .lt. 0.d0) then	! use # contours
	    cd = (fmin-fmax)/ev.dfdoa
	  endif
	  cs = INT(fmin/cd)*cd
	  nn = (fmax-cs)/cd + 1
	elseif (ev.fdmode .eq. PG_REGION) then
	  cs = -ev.dfdoa/2
	  cd = ev.dfdoa
	  nn = 2
	elseif (ev.fdmode .eq. PG_LINE) then
	  cs = 0
	  cd = ev.dfdoa
	  nn = 1
	else
	  nn = 0
	endif
	if (nn.gt.0) then
	  call MX$CONTOUR (cs,cd,nn,fd,pchk,
     &      buf(7), lnf, L_ctmp, mg,mg,cxs,cxd,cys,cyd)
	  if (lnf.gt.0) then
	    ly.nelem = lnf+3
	    ly.color = PC_FDOA
	    call XGEO_DRAW_LAYER (ly)
	    Gx.evp(5).lon = buf(13) + xs
	    Gx.evp(5).lat = buf(14) + ys
	  endif
	endif
c
c draw the horizon lines
	if (BTEST(ev.flags,PM_HORIZON)) then
	  ly.color = PC_HORIZON
	  if (Gx.svh) ly.color = ly.color+(ev.iballa-1)
	  call MX$CONTOUR (0.0,0.0,1,aa,pchk,
     &      buf(7), lne, L_ctmp, mg,mg,cxs,cxd,cys,cyd)
	  ly.nelem = lne+3
	  call XGEO_DRAW_LAYER (ly)
	  call MX$CONTOUR (0.0,0.0,1,ab,pchk,
     &      buf(7), lne, L_ctmp, mg,mg,cxs,cxd,cys,cyd)
	  ly.nelem = lne+3
	  call XGEO_DRAW_LAYER (ly)
	endif
c
c restore line width
	call MX$DRAW_MODE (1)
c
c draw points
	if (BTEST(ev.flags,PM_POINTS)) then
	  ly.type = 5000
	  ly.name = '*'
	  stmp = CHAR( n+48 )
	  ll = 2
	  if (Gx.nevent.le.1) ll = 1
	  Gx.evp(1).name = 'S'//stmp
	  Gx.evp(1).name_len = ll
	  call SET0 (ev.site,Gx.evp(1))
	  if (ev.svmode .eq. PG_CURRENT) then
	    Gx.evp(2).name = 'V'//CHAR(ev.iballa+48)
	    Gx.evp(3).name = 'V'//CHAR(ev.iballb+48)
	  else
	    Gx.evp(2).name = 'A'//stmp
	    Gx.evp(2).name_len = ll
	    Gx.evp(3).name = 'B'//stmp
	    Gx.evp(3).name_len = ll
	  endif
	  call SET0 (ev.sva,Gx.evp(2))
	  call SET0 (ev.svb,Gx.evp(3))
	  Gx.evp(4).name = 'T'//stmp
	  Gx.evp(4).name_len = ll
	  Gx.evp(5).name = 'F'//stmp
	  Gx.evp(5).name_len = ll
	  ly.nelem = 5
	  if (ev.fdmode .eq. PG_NONE) then
	    ly.nelem = ly.nelem-1
	  endif
	  if (ev.tdmode .eq. PG_NONE) then
	    Gx.evp(4) = Gx.evp(5)
	    ly.nelem = ly.nelem-1
	  endif
	  do i = 1,ly.nelem
	    Gx.evp(i).alt = 0
	  enddo
	  ly.color = -1
	  ly.bufloc = %loc(Gx.evp)
	  call XGEO_DRAW_LAYER (ly)
	endif

	return
	end

	subroutine SET0 (s,t)
	real*8 s(3),t(3),x,y,z,r,f1,a
	include 'constants.inc'
	f1  = 1.D0 - earth_flattening_factor
	f2  = (2.D0-earth_flattening_factor)*earth_flattening_factor 
	x = s(1)
	y = s(2)
	z = s(3)
	r = sqrt (x*x+y*y+z*z)
	a = asin(z/r)
	ca = cos(a)
	t(1) = r - (earth_equatorial_radius*f1)/dsqrt(1.d0-f2*ca*ca)
	t(2) = drad2deg * atan( tan(a)/(f1*f1) )
	t(3) = drad2deg * atan2(y,x)
	return
	end

	subroutine GEO2ECR (alt,lat,lon,xt,yt,zt,rr)
	implicit real*8 (A-Z)
	include 'constants.inc'
c
        ! R = ERADEQ * F1/SQRT(1-F2*CPHI*CPHI)
	f1  = 1.D0 - earth_flattening_factor
	f2  = (2.D0-earth_flattening_factor)*earth_flattening_factor 
c
	rlat = lat*ddeg2rad
	rlat = datan(dtan(rlat)*f1*f1)
	rlon = lon*ddeg2rad
	slat = dsin(rlat)
	clat = dcos(rlat)
	slon = dsin(rlon)
	clon = dcos(rlon)
	rr = (earth_equatorial_radius*f1)/dsqrt(1.d0-f2*clat*clat) + alt
	xt = rr*clat*clon
	yt = rr*clat*slon
	zt = rr*slat
c
	return
	end



	subroutine ECR2GEO (xt,yt,zt,lon,lat,alt)
	implicit real*8 (A-Z)
	include 'constants.inc'
c
	call M$ERROR ('Not coded yet')
c
        ! R = ERADEQ * F1/SQRT(1-F2*CPHI*CPHI)
	f1  = 1.D0 - earth_flattening_factor
	f2  = (2.D0-earth_flattening_factor)*earth_flattening_factor 
c
	rlat = lat*ddeg2rad
	rlat = datan(dtan(rlat)*f1*f1)
	rlon = lon*ddeg2rad
	slat = dsin(rlat)
	clat = dcos(rlat)
	slon = dsin(rlon)
	clon = dcos(rlon)
	rr = (earth_equatorial_radius*f1)/dsqrt(1.d0-f2*clat*clat) + alt
	xt = rr*clat*clon
	yt = rr*clat*slon
	zt = rr*slat
c
	return
	end



	subroutine TRIM_NAME (filename, target)
c
	character*(*)   filename,target
	integer*4       i,j
c
c trim name of path and extensions
	i = M$SEARCHB(filename,']')
	if (i.eq.0) i = M$SEARCHB(filename,'/')
	if (i.eq.0) i = M$SEARCHB(filename,':')
	j = M$SEARCH(filename(i+1:),'.') - 1
	if (j.lt.0) j = M$LENGTH(filename) - i
	target = filename(i+1:i+j)
	call M$UPPERCASE (target(1:j))
c
	return
	end

	subroutine XGEO_GET_SVS (ev,time)
	include 'xgeoplot.inc'
	if (Gx.nevent .eq. 0) return
	ev.time = time
	if (ev.gsmode.ne.PG_FIXED) call M$MOVE (Gx.site, ev.site, 24)
	if (Gx.nball .eq. 0) return
	if (ev.iballa .le. Gx.nball .and. ev.iballb .le. Gx.nball) then
	  call XGEO_GET_SV (%val(Gx.lhsv(ev.iballa)), 
     &		ev.sva, ev.time, ev.site, BTEST(ev.flags,PM_LTC))
	  call XGEO_GET_SV (%val(Gx.lhsv(ev.iballb)), 
     &		ev.svb, ev.time, ev.site, BTEST(ev.flags,PM_LTC))
	else
	  call M$BEEP (1,0.0)
	  call M$WARNING ('Event specifies a SV for ball > NBALL')
	endif
	return
	end


	subroutine XGEO_GET_SV (hsv, sv, time, site, ltc)
	include 'headers.inc'
	include 'constants.inc'
	record /HEADER/ hsv
	integer*4 nget
	logical*4 ltc
	real*8 dval,time,sv(6),site(3),svt(6,4),ts,tt

	if (time.lt.86400*100) then	! no date in time
	  tt = time
	elseif (hsv.timecode .ne. 0) then 
	  tt = time - hsv.timecode	! use time ref to header start
	else
	  tt = mod (time, 86400.d0)	! use time modulo one day
	endif

	dval = DNINT( (tt - hsv.tstart)/hsv.tdelta) + 1
	dval = max(1.d0, min(hsv.size-3.d0, dval-1) )
	nget = min(4,INT(hsv.size))
	call M$GRAB5N (hsv,svt,dval,nget,n,'POS,VEL')
	if (n.ne.nget) call MX$MESSAGE ('Bad state vector read',1.0)

	ts = hsv.tstart + hsv.tdelta*(dval-1)
	call M$INTERPD (nget,ts,hsv.tdelta,svt, 1,tt,0.d0,sv, 6,3,n)

	if (ltc) then
	  dval = dsqrt( (sv(1)-site(1))**2 + (sv(2)-site(2))**2 +
     &		(sv(3)-site(3))**2 ) / speed_of_light
	  call M$INTERPD (nget,ts,hsv.tdelta,svt, 1,tt-dval,0.d0,sv, 6,3,n)
	endif

	return
	end


	subroutine XGEO_ALTITUDE (hcb, p)
	include 'xgeoplot.inc'
	record /HEADER/ hcb
	record /GEOPOS/ p
	real*8 s
	ix = INT((p.x-hcb.xstart)/hcb.xdelta)
	iy = INT((p.y-hcb.ystart)/hcb.ydelta)
	s = iy*hcb.subsize + ix + 1
	if (s.ge.1 .and. s.le.hcb.size .and. 
     +	    ix.gt.0 .and. ix.lt.hcb.subsize) then
	  call M$GRAB (hcb,p.z,s,1,-1)
	else
	  p.z = -1
	endif
	return
	end


	subroutine PREDICT_EVAL (m,ev, alt,lat,lon, tdoa,fdoa, aoaa,aoab)
c
	include 'xgeoplot.inc'
	real*8 lat,lon,alt,tdoa(3),fdoa(3),aoaa,aoab, x,y,z
	real*8 xt,yt,zt,rr,rri,cci,scale
	real*8 xta,yta,zta,rta,rtai,rrta,rrtar
	real*8 xtb,ytb,ztb,rtb,rtbi,rrtb,rrtbr
	real*8 xsa,ysa,zsa,rsa,rrsa
	real*8 xsb,ysb,zsb,rsb,rrsb
	real*8 f1,f2,rlat,rlon,slat,slon,clat,clon
	real*8 dxdlat,dydlat,dzdlat,dxdlon,dydlon,dzdlon
	real*8 dtddx,dtddy,dtddz,dfddx,dfddy,dfddz
	real*8 dtddlat,dfddlat,dtddlon,dfddlon
c
c convert geodetic lat,lon,alt to ECR cartesian
c	call GEO2ECR (alt,lat,lon,xt,yt,zt,rr)
c
        ! R = ERADEQ * F1/SQRT(1-F2*CPHI*CPHI)
	f1  = 1.D0 - earth_flattening_factor
	f2  = (2.D0-earth_flattening_factor)*earth_flattening_factor 
	rlat = lat*ddeg2rad
	rlat = datan(dtan(rlat)*f1*f1)
	rlon = lon*ddeg2rad
	slat = dsin(rlat)
	clat = dcos(rlat)
	slon = dsin(rlon)
	clon = dcos(rlon)
	rr = (earth_equatorial_radius*f1)/dsqrt(1.d0-f2*clat*clat) + alt
	rri = 1.d0/rr
	cci = 1.d0/ speed_of_light
	xt = rr*clat*clon
	yt = rr*clat*slon
	zt = rr*slat
	goto 2

	entry PREDICT_EVAL_ECR (m,ev, x,y,z, tdoa,fdoa, aoaa,aoab)
	xt = x
	yt = y
	zt = z
	rri  = 1.0d0/dsqrt(xt*xt + yt*yt + zt*zt)

    2	continue
c
c calculate relative distances to receivers
	xta  = ev.sva(1) - xt
	yta  = ev.sva(2) - yt
	zta  = ev.sva(3) - zt
	xtb  = ev.svb(1) - xt
	ytb  = ev.svb(2) - yt
	ztb  = ev.svb(3) - zt
c
c calculate range to receivers
	rta  = dsqrt ( xta*xta + yta*yta + zta*zta )		
	rtb  = dsqrt ( xtb*xtb + ytb*ytb + ztb*ztb )		
	rtai = 1.d0/rta
	rtbi = 1.d0/rtb
c
c calculate sines of angle of arrival values
	aoaa = (xt*xta + yt*yta + zt*zta) * rtai * rri
	aoab = (xt*xtb + yt*ytb + zt*ztb) * rtbi * rri
c
c return if only interested in AOA
	if (m.lt.0) return
c
c calculate range rates to receivers
	rrta = ( xta*ev.sva(4) + yta*ev.sva(5) + zta*ev.sva(6) ) * rtai		
	rrtb = ( xtb*ev.svb(4) + ytb*ev.svb(5) + ztb*ev.svb(6) ) * rtbi		
c
c calculate relative distances from site to receivers
	xsa  = ev.sva(1) - ev.site(1)
	ysa  = ev.sva(2) - ev.site(2)
	zsa  = ev.sva(3) - ev.site(3)
	xsb  = ev.svb(1) - ev.site(1)
	ysb  = ev.svb(2) - ev.site(2)
	zsb  = ev.svb(3) - ev.site(3)
c
c calculate range and rangerates from site to receivers at t = time
	if (BTEST(ev.flags,PM_DWNLNK)) then 
	  rsa  = dsqrt ( xsa*xsa + ysa*ysa + zsa*zsa )		
	  rsb  = dsqrt ( xsb*xsb + ysb*ysb + zsb*zsb )		
	  rrsa = ( xsa*ev.sva(4) + ysa*ev.sva(5) + zsa*ev.sva(6) ) / rsa
	  rrsb = ( xsb*ev.svb(4) + ysb*ev.svb(5) + zsb*ev.svb(6) ) / rsb
	else		! no downlink
	  rsa = 0
	  rsb = 0
	  rrsa = 0
	  rrsb = 0
	endif
c
c calculate TDOA and FDOA
	tdoa(1) = cci* ( rtb - rta + rsb - rsa )	
	fdoa(1) = -cci* ( ev.flinkt*(rrtb-rrta) + 
     &		ev.flinkb*rrsb - ev.flinka*rrsa )	
c
c return if not interested in 1st derivative measurements
	if (m.eq.0) return
	rrtar = rrta*rtai
	rrtbr = rrtb*rtbi
c
c calculate derivatives x,y,z vs. lat,lon
	dxdlat = -slat*clon
	dydlat = -slat*slon
	dzdlat = clat
	dxdlon = -clat*slon
	dydlon = clat*clon
	dzdlon = 0
c
c calculate derivatives tdoa,fdoa vs. x,y,x 
	dtddx = (xta*rtai) - (xtb*rtbi)
	dtddy = (yta*rtai) - (ytb*rtbi)
	dtddz = (zta*rtai) - (ztb*rtbi)
	dfddx = (xta*rrtar-ev.sva(4))*rtai - (xtb*rrtbr-ev.svb(4))*rtbi
	dfddy = (yta*rrtar-ev.sva(5))*rtai - (ytb*rrtbr-ev.svb(5))*rtbi
	dfddz = (zta*rrtar-ev.sva(6))*rtai - (ztb*rrtbr-ev.svb(6))*rtbi
c
c calculate derivatives tdoa,fdoa vs. lat,lon
	dtddlat = dtddx*dxdlat + dtddy*dydlat + dtddz*dzdlat
	dtddlon = dtddx*dxdlon + dtddy*dydlon + dtddz*dzdlon
	dfddlat = dfddx*dxdlat + dfddy*dydlat + dfddz*dzdlat
	dfddlon = dfddx*dxdlon + dfddy*dydlon + dfddz*dzdlon
	scale = rr*cci*ddeg2rad  ! in units sec/deg
	tdoa(2) = dtddlat*scale
	tdoa(3) = dtddlon*scale
	scale = scale*ev.flinkt  ! in units Hz/deg
	fdoa(2) = dfddlat*scale
	fdoa(3) = dfddlon*scale
c
c return if derivatives w.r.t. degrees
	if (m.eq.1) return
c
c calculate derivatives w.r.t. meters
	scale = 1.d0 / (ddeg2rad*rr)
	tdoa(2) = tdoa(2) * scale
	fdoa(2) = fdoa(2) * scale
	scale = scale / clat 
	tdoa(3) = tdoa(3) * scale
	fdoa(3) = fdoa(3) * scale
c
	return
	end


	subroutine XGEO_SET_PERSP (iball)
	include 'xgeoplot.inc'
	integer*4 iball
	real*8 sv(6)

	call XGEO_GET_SV (%val(Gx.lhsv(iball)), 
     &			sv, Gx.time, Gx.site, 0)
	dval = dsqrt( sv(1)**2 + sv(2)**2 + sv(3)**2 )
	Gx.plat = drad2deg*dasin(sv(3)/dval)
	Gx.plon = drad2deg*datan2(sv(2),sv(1))
	Gx.pdis = dval/earth_equatorial_radius

	return
	end


	subroutine XGEO_COLORS (mm,str)
	include 'xgeoplot.inc'
	record /HEADER/ hcb

	integer*2 map(4,L_colors) /
     &	-1,30,30,30, -1,0,80,0,   -1,80,0,0,   -1,80,80,0,  
     &	-1,80,0,80,  -1,30,40,55, -1,45,45,45, -1,35,30,30, 
     &	-1,0,20,50,  -1,20,25,20, -1,40,40,0,  -1,99,99,90, 
     &	-1,100,50,0, -1,100,80,70,-1,0,50,100, -1,100,100,100,
     &	-1,0,0,40,   -1,0,6,42,  -1,0,13,44,  -1,0,19,46,
     &	-1,0,26,48,  -1,0,32,50,  -1,0,38,52,  -1,0,45,54,
     &	-1,34,81,45, -1,43,72,41, -1,51,63,36, -1,59,55,32,
     &	-1,68,46,28, -1,78,52,38, -1,88,74,64, -1,100,100,95  /

	if (mm.eq.0) then
	  if (str.eq.' ') then
	    do i = 1,Gx.ncolors
	      call MX$COLORTAB (i,map(1,i),-2)
	    enddo
	  else
 	    call M$INIT (hcb,str,'1000','QI',0)
	    call M$OPEN (hcb,HCBF_INPUT+HCBF_NOABORT)
	    if (hcb.open) then
	      n = min(NINT(hcb.size),Gx.ncolors)
	      call M$GRAB (hcb,map,1.d0,n,-1)
	      call M$CLOSE (hcb)
	      do i = 1,n
		call MX$COLORTAB (i,map(1,i),-2)
	      enddo
	    else
	      call MX$MESSAGE ('Error opening colormap file',-1.0)
	    endif
	  endif

	elseif (mm.eq.-1) then

 2	  m = MX$MENU (1,'COLORMAP;OK,Alter,Load,Save')
	  if (m.eq.2) then
	    n = MX$MENU (0,colors)
	    if (n.eq.12) n = MX$MENU(0,'USERMAP;1,2,3,4,5,6,7,8,9,10,'//
     +		'11,12,13,14,15,16,17,18,19,20') + 12
	    if (n.eq.0 .or. n.eq.12) goto 2
	    Mx.xpos = Mx.width/2
	    Mx.ypos = Mx.height/2
	    call MX$COLORTAB (n,map,-1)
	    call MX$COLORTAB (n,map(1,n),-3)
	  elseif (m.eq.3) then
	    hcb.file_name = ' '
	    call MX$PROMPT ('Filename: ',hcb.file_name)
 	    call M$INIT (hcb,hcb.file_name,'1000','QI',0)
	    call M$OPEN (hcb,HCBF_INPUT+HCBF_NOABORT)
	    if (hcb.open) then
	      n = min(NINT(hcb.size),Gx.ncolors)
	      call M$GRAB (hcb,map,1.d0,n,-1)
	      call M$CLOSE (hcb)
	      do i = 1,n
		call MX$COLORTAB (i,map(1,i),-2)
	      enddo
	    else
	      call MX$MESSAGE ('Error opening colormap file',-1.0)
	    endif
	  elseif (m.eq.4) then
	    hcb.file_name = ' '
	    call MX$PROMPT ('Filename: ',hcb.file_name)
	    call M$INIT (hcb,hcb.file_name,'1000','QI',0)
	    hcb.size = Gx.ncolors
	    call M$OPEN (hcb,HCBF_OUTPUT+HCBF_NOABORT)
	    if (hcb.open) then
	      call M$FILAD (hcb,map,Gx.ncolors)
	      call M$CLOSE (hcb)
	    else
	      call MX$MESSAGE ('Error opening colormap file',-1.0)
	    endif
	  endif
	  if (m.ne.1) goto 2
	endif
	return
	end



	subroutine XGEO_EVENTS_SAVE (name,k)
	include 'xgeoplot.inc'
	record /HEADER/ hcb
	character*(*) name
	pointer (pev,ev)
c
	if (name .eq. ' ') return
	call M$INIT (hcb,name,'5011','NH',0)
	n = M$PUT_KEYDATA (hcb,'INDEX',k,4,0,'L')
	hcb.frame_of_ref = 'ECR'
	call M$ADDCOMP (hcb,'NAME','4A',0,0)
	call M$ADDCOMP (hcb,'MODE','8B',0,0)
	call M$ADDCOMP (hcb,'FLAG','2L',0,0)
	call M$ADDCOMP (hcb,'TIME','2D',1,1)
	call M$ADDCOMP (hcb,'TDOA','2D',1,1)
	call M$ADDCOMP (hcb,'FDOA','2D',1,3)
	call M$ADDCOMP (hcb,'FRQT','SD',1,3)
	call M$ADDCOMP (hcb,'FRQA','SD',1,3)
	call M$ADDCOMP (hcb,'FRQB','SD',1,3)
	call M$ADDCOMP (hcb,'POSS','VD',2,5)
	call M$ADDCOMP (hcb,'POSA','VD',2,5)
	call M$ADDCOMP (hcb,'VELA','VD',2,6)
	call M$ADDCOMP (hcb,'POSB','VD',2,5)
	call M$ADDCOMP (hcb,'VELB','VD',2,6)
	hcb.size = Gx.nevent
	call M$OPEN (hcb,HCBF_OUTPUT+HCBF_NOABORT)
	do n = 1,Gx.nevent
	  pev = Gx.eloc(n)
	  call M$FILAD (hcb, ev, 1)
	enddo
	call M$CLOSE (hcb)
	return

	entry XGEO_EVENTS_LOAD (name,k)

	if (name .eq. ' ') return
	call M$INIT (hcb,name,'5011','NH',0)
	call M$OPEN (hcb,HCBF_INPUT+HCBF_NOABORT)
	if (.not. hcb.open) then
	  call MX$MESSAGE ('Error opening file: '//hcb.file_name,-1.0)
	  return
	endif
	if (hcb.size .gt. L_eventmx-Gx.nevent) then
	  call MX$MESSAGE ('Too many events - truncating',1.0)
	  hcb.size = L_eventmx-Gx.nevent
	endif
	if (k.eq.-1 .and. hcb.size.gt.0) call XGEO_UNLOAD (P_EVENT)
	do n = 1,hcb.size
	  Gx.nevent = Gx.nevent + 1
	  call M$MALLOC (L_eventsz, Gx.eloc(Gx.nevent) )
	  call XGEO_EVENT_INIT (Gx.nevent)
	  pev = Gx.eloc(Gx.nevent)
	  call M$GRAB (hcb, ev, DBLE(n), 1, -1)
	enddo
	call M$CLOSE (hcb)
	return

	end


	function XGEO_LOCATE ()
	include 'xgeoplot.inc'
	record /EVENT/ ev1,ev2
	pointer (pev1,ev1),(pev2,ev2)
c
	real*8 da(3),ta(3),fa(3)
	real*8 aoaa,aoab,tolr,tdoa,fdoa
	real*8 xa,ya,za,dx,dy,dt,df,dtx,dty,dfx,dfy,dlmin,dlmax
	character*2 ss
c
	pev1 = Gx.eloc(Gx.glevn1)
	pev2 = Gx.eloc(Gx.glevn2)
c
	dlmin = -9.0
	dlmax = 11.0
	tolr = 1.e-8
	ii = 0
c
  2	ii = ii+1
	xa = Gx.mrk.x 
	ya = Gx.mrk.y 
	za = Gx.mrk.z 
c
	if (ii.gt.100) then
	  XGEO_LOCATE = 0
	  return
	endif
c
	if (Gx.glmode.eq.1 .and. Gx.glevn1.eq.Gx.glevn2) then
	  call PREDICT_EVAL (1,ev1, za,ya,xa, ta,fa,aoaa,aoab)
	  tdoa = ev1.tdoa
	  fdoa = ev1.fdoa
	elseif (Gx.glmode.eq.1) then
	  call PREDICT_EVAL (1,ev1, za,ya,xa, ta,da,aoaa,aoab)
	  tdoa = ev1.tdoa
	  call PREDICT_EVAL (1,ev2, za,ya,xa, da,fa,aoaa,aoab)
	  fdoa = ev2.fdoa
	elseif (Gx.glmode.eq.2) then
	  call PREDICT_EVAL (1,ev1, za,ya,xa, ta,da,aoaa,aoab)
	  tdoa = ev1.tdoa
	  call PREDICT_EVAL (1,ev2, za,ya,xa, fa,da,aoaa,aoab)
	  fdoa = ev2.tdoa
	elseif (Gx.glmode.eq.3) then
	  call PREDICT_EVAL (1,ev1, za,ya,xa, da,ta,aoaa,aoab)
	  tdoa = ev1.fdoa
	  call PREDICT_EVAL (1,ev2, za,ya,xa, da,fa,aoaa,aoab)
	  fdoa = ev2.fdoa
	endif
c
	dt = tdoa - ta(1)
	df = fdoa - fa(1)
	dtx = ta(3)
	dty = ta(2)
	dfx = fa(3)
	dfy = fa(2)
	dy = (dt*dfx-df*dtx)/(dfx*dty-dtx*dfy)
	dx = (df-dfy*dy)/dfx
c
	if (BTEST(Mc.debug,5)) then
	   write(*,*) 'locate ',ii,dx,dy
	   k = Mx.level
	   ix = (Gx.mrk.x-Mx.stk(k).xmin)/Mx.stk(k).xscl + Mx.stk(k).x1
	   iy = (Mx.stk(k).ymax-Gx.mrk.y)/Mx.stk(k).yscl + Mx.stk(k).y1
	   k = M$L2A (ii,ss)
	   call MX$TEXT (ix,iy,ss)
	endif
c
	if (dy.gt.dlmax) then
	  dy = dlmax
	elseif (dy.lt.dlmin) then
	  dy = dlmin
	endif
	Gx.mrk.y = Gx.mrk.y + dy
	if (Gx.mrk.y .lt. -88.) then
	  Gx.mrk.y = max(-88.d0, -180.d0-Gx.mrk.y)
	  Gx.mrk.x = Gx.mrk.x+180.d0 
	elseif (Gx.mrk.y .gt. 88.) then
	  Gx.mrk.y = min(88.d0, 180.d0-Gx.mrk.y)
	  Gx.mrk.x = Gx.mrk.x+180.d0 
	endif
c
	if (dx.gt.dlmax) then
	  dx = dlmax
	elseif (dx.lt.dlmin) then
	  dx = dlmin
	endif
	Gx.mrk.x = Gx.mrk.x + dx
	if (Gx.mrk.x .lt. -180.+Gx.xoffs) then
	  Gx.mrk.x = Gx.mrk.x+360.d0 
	elseif (Gx.mrk.x .gt.  180.+Gx.xoffs) then
	  Gx.mrk.x = Gx.mrk.x-360.d0 
	endif
c
	if (dabs(dx).gt.tolr) goto 2
	if (dabs(dy).gt.tolr) goto 2
c
	XGEO_LOCATE = 1
	if (BTEST(Mc.debug,5)) i = mx$menu(1,'ok')
c
	return
	end


	subroutine CEP_EVAL (alt,lat,lon, cep)
	include 'xgeoplot.inc'
	record /EVENT/ ev1,ev2
	pointer (pev1,ev1),(pev2,ev2)
	real*8 alt,lat,lon,cep,ta(3),fa(3),da(3),aoaa,aoab,ellipse(3)

	pev1 = Gx.eloc(Gx.glevn1)
	pev2 = Gx.eloc(Gx.glevn2)

	if (Gx.glmode.eq.1 .and. pev1.eq.pev2) then
	  call PREDICT_EVAL (2,ev1,alt,lat,lon, ta,fa, aoaa,aoab)
	elseif (Gx.glmode.eq.1) then
	  call PREDICT_EVAL (2,ev1,alt,lat,lon, ta,da, aoaa,aoab)
	  call PREDICT_EVAL (2,ev2,alt,lat,lon, da,fa, aoaa,aoab)
	elseif (Gx.glmode.eq.2) then
	  call PREDICT_EVAL (2,ev1,alt,lat,lon, ta,da, aoaa,aoab)
	  call PREDICT_EVAL (2,ev2,alt,lat,lon, fa,da, aoaa,aoab)
	elseif (Gx.glmode.eq.3) then
	  call PREDICT_EVAL (2,ev1,alt,lat,lon, da,ta, aoaa,aoab)
	  call PREDICT_EVAL (2,ev2,alt,lat,lon, da,fa, aoaa,aoab)
	endif
	call ELLIPSE_CEP (ta,fa, Gx.glstd1,Gx.glstd2, ellipse,cep)
	return
	end

	subroutine ELLIPSE_CEP (tvec,fvec, TDOA_std,FDOA_std, ellipse,cep)
c
c
c	Tvec(1)		contains the TDOA value in seconds (not used)
c	Tvec(2)		contains the TDOA partial w.r.t. North in sec/meter
c	Tvec(3)		contains the TDOA partial w.r.t. East in sec/meter
c
c	Fvec(1)		contains the FDOA value in Hz (not used)
c	Fvec(2)		contains the FDOA partial w.r.t. North in Hz/meter
c	Fvec(3)		contains the FDOA partial w.r.t. East in Hz/meter
c
c	TDOA_std	contains the Standard Deviation of the TDOA in Seconds
c	FDOA_std	contains the Standard Deviation of the FDOA in Hertz
c
c	Output:  All Estimates are 95% confidence
c
c	ELLIPSE(1)	contains the Semi-Major Axis of the Ellipse (meters)
c	ELLIPSE(2)	contains the Semi-Minor Axis of the Ellipse (meters)
c	ELLIPSE(3)	contains the Major Axis of the Ellipse (degrees)
c
c	CEP		contains the circular error prob. (meters)
c
	real*8 Tvec(3), Fvec(3)
	real*4 TDOA_std, FDOA_std
	real*8 ellipse(3),cep
c
c
	real*8 sensitive(2,2)
	real*8 covariance(2,2)
	real*8 TDOA_var, FDOA_var
	real*8 det, deti
	real*8 temp1,temp2,tp1,tp2
	real*8 m, v, ratio

c
c	   From Geometry; (J from R. Brubaker's paper)
c
c		|			|
c		| T(North)	T(East)	|
c	J   =	|			|
c		| F(North)	F(East)	|
c		|			|
c
	det = Tvec(2)*Fvec(3) - Tvec(3)*Fvec(2)
	det = SIGN( max(1d-22,abs(det)), det )
c	deti = 1.0d0/det
	deti = 1.0d-6/det	 	! to prevent numerical overflow
c
c	Calculate the inverse 
c	(The sensistivity Matrix is called J inverse in Brubaker's paper)
c
	sensitive(1,1) =  Fvec(3)*deti
	sensitive(2,2) =  Tvec(2)*deti
	sensitive(1,2) = -Tvec(3)*deti
	sensitive(2,1) = -Fvec(2)*deti

	TDOA_var = TDOA_std * TDOA_std 
	FDOA_var = FDOA_std * FDOA_std 
c
c	Covariance Matrix is:  sensitive * T/C * sensitive transpose
c
	covariance(1,1) = TDOA_var * sensitive(1,1) * sensitive(1,1) +
     &			  FDOA_var * sensitive(1,2) * sensitive(1,2)

	covariance(2,2) = TDOA_var * sensitive(2,1) * sensitive(2,1) +
     &			  FDOA_var * sensitive(2,2) * sensitive(2,2) 

	covariance(1,2) = TDOA_var * sensitive(1,1) * sensitive(2,1) +
     &			  FDOA_var * sensitive(1,2) * sensitive(2,2) 

	covariance(2,1) = covariance(1,2) 

	temp1 = (covariance(1,1) +  covariance(2,2)) * 0.5
	temp2 = sqrt( ( covariance(1,1)-covariance(2,2) )**2 * .25
     &		+  covariance(1,2)**2 )

	ellipse(1) = 2.4477 * sqrt( temp1 + temp2 )
	ellipse(2) = 2.4477 * sqrt( temp1 - temp2 )
	ellipse(3) = .5* atan2d( 2.0 * covariance(1,2), 
     &			covariance(1,1) - covariance(2,2) )
	if( ellipse(3) .lt. 0 ) ellipse(3) = ellipse(3) + 180.0d0

	tp1 = temp1 + temp2
	tp2 = temp1 - temp2
	m = tp1 + tp2
	v = 2.0 * ( tp1*tp1 + tp2*tp2 )
	ratio = v / (9 * m * m)
	
	cep = sqrt( m * ( 1.6447d0 * sqrt(ratio) 
     &		        + (1.0d0 - ratio) )**3.0d0 )
c	cep = cep * 1.0d0 /1.852d3	! in nautical miles
c	cep = cep * 1.0d6 /1.852d3	! to reverse prevent numerical overflow
	cep = cep * 1.0d6 ! to reverse prevent numerical overflow
	return
	end

	subroutine XGEO_DUMP_PREDICTS (ev,t1arg,t2arg,td,name)
	include 'xgeoplot.inc'
	record /HEADER/ hcb,hsv
	real*8 td,t1arg,t2arg,rec(15)
	character*(*) name
	pointer (phsv,hsv)
c
	phsv = Gx.lhsv(ev.iballa) 
	if (t1arg.lt.86400*100) then	! no date in time
	  t1 = t1arg
	  t2 = t2arg
	elseif (hsv.timecode .ne. 0) then 
	  t1 = t1arg - hsv.timecode	! use time ref to header start
	  t2 = t2arg - hsv.timecode	! use time ref to header start
	else
	  t1 = mod (t1arg, 86400.d0)	! use time modulo one day
	  t2 = mod (t2arg, 86400.d0)	! use time modulo one day
	endif
c
	call M$INIT (hcb,name,'5000','NH',0)
	hcb.timecode = hsv.timecode
	hcb.tstart = t1
	hcb.tdelta = td
	hcb.frame_of_ref = 'TOPOCENT'
	hcb.altitude  = Gx.mrk.alt
	hcb.latitude  = Gx.mrk.lat
	hcb.longitude = Gx.mrk.lon
	call M$ADDCOMP (hcb,'POS ','VD',2,5)
	call M$ADDCOMP (hcb,'TDOA','SD',1,1)
	call M$ADDCOMP (hcb,'DTDY','SD',1,1)
	call M$ADDCOMP (hcb,'DTDX','SD',1,1)
	call M$ADDCOMP (hcb,'FDOA','SD',1,3)
	call M$ADDCOMP (hcb,'DFDY','SD',1,3)
	call M$ADDCOMP (hcb,'DFDX','SD',1,3)
	call M$ADDCOMP (hcb,'ELPS','VD',0,5)
	call M$ADDCOMP (hcb,'CEP ','SD',1,5)
	call M$ADDCOMP (hcb,'AOAA','SD',0,34)
	call M$ADDCOMP (hcb,'AOAB','SD',0,34)
c
	rec(1) = 0
	rec(2) = 0
	rec(3) = 0
	hcb.size = 0
	if (td.gt.0.d0) hcb.size = NINT((t2-t1)/td) + 1
	call M$OPEN (hcb,HCBF_OUTPUT+HCBF_NOABORT)
	do time = t1,t2,td
	  call XGEO_GET_SVS (ev,time)
	  call PREDICT_EVAL (2,ev,Gx.mrk.alt,Gx.mrk.lat,Gx.mrk.lon,
     &		rec(4),rec(7),rec(14),rec(15))
	  call ELLIPSE_CEP (rec(4),rec(7), Gx.glstd1,Gx.glstd2, 
     &		rec(10),rec(13))
	  call M$FILAD (hcb,rec,1)
	enddo
	call XGEO_GET_SVS (ev,Gx.time)
c
	call M$CLOSE (hcb)
	return
	end

	subroutine T5INIT (t5, mode, color, symbol, size)
	include 'xgeoplot.inc'
	record /T5REC/ t5
	integer*4 mode,color,symbol,size
	t5.x = 0
	t5.y = 0
	t5.z = 0
	t5.name_len = 0
	t5.name_pos = 0
	t5.color  = color
	t5.colorbg= 0
	t5.cone   = mode
	t5.invert = 0
	t5.symbol = symbol
	t5.size   = size
	return
	end


	subroutine M$REMALLOC_EXT (size, location, needsize)
	integer*4 size, location, needsize, newlocation, newsize
	if (needsize.gt.size) then
	  newsize = needsize + 1024
	  call M$MALLOC (newsize, newlocation)
	  call M$MOVE (%val(location), %val(newlocation), size)
	  call M$FREE (location)
	  location = newlocation
	  size = newsize
	endif
	return
	end
