+title.
      SPK      2.00 /00      990819   2136

   ***********************************************************
   **  Spark Chamber Simulation and Reconstruction Program  **
   ***********************************************************

   History

   29.07.98   FC   1.00/00 created from ZUT graphics for spark chamber
   19.08.99   FC   2.00/00 new data without reconstruction, annoying bugs
   19.08.99   FC   distance curves put in, work right away, suspicious
   20.08.99   FC   also Monte-Carlo data in, force known points in
   20.08.99   FC   call igfais(3) -> igset('fais',0) solves it, why?
   23.08.99   FC   add the spkcir routine, average per level for data
   24.08.99   FC   generalize mc and reconstruction parts
   25.08.99   FC   last step with trajectories if enough info (e.g. MC)
   25.08.99   FC   cleanup, comment, put in minuscules
   01.08.02   FC   put plot arguments in upper, e.g. 'L', 'S', ..
   ++++++++

+patch,spkcom.   The SPK common blocks
+keep,spkios.
c  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c
c  general and i/o global declarations
c
c     constants (              rad,       rad/deg)
      parameter (pi=3.141592653587,pi180=pi/180.0)
c
      common/spkios/       ldate ,ltime ,
     +              lvers ,lverd ,lvert ,
     +              lunin ,lunout,lunhpl,messag
c
c  where:  ldate  = current date in the yymmdd integer format
c          ltime  = .       time .  .   hhmm   .       .
c          lvers  = version number e.g. 10504
c          lverd  = .       date in the yymmdd integer format
c          lvert  = .       time .  .   hhmm   .       .
c
c          lunin  = input  unit number for changes to initialization
c          lunout = output .    .      .   messages, printout, etc..
c          lunhpl = output .    .      .   plots
c                                                               routines
c          messag = output message level (-1:none      )
c
+keep,spkgeo.
c  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c
c  spark chamber geometry
c
      common/spkgeo/platex,platey,spacez,thickz,nplate,bounds(3)
c
c          platex = half-size of      plates in x [cm]  (horizontal)
c          platey = half-size of      plates in y [cm]  (horizontal)
c          spacez = spacing between   plates in z [cm]  (vertical)
c          thickz = half-thickness of plates in z [cm]
c          nplate = number of plates (central one at z=0)
c          bounds = half-limits of standard displays (x/y/z) [cm]
c
+keep,spkdat.
c  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c
c  spark chamber data
c
      parameter (mmic=8,mlev=2)
      common/spkdat/nmic,posi(3,mmic),data(mmic),dist(mmic),lepo(mmic),
     +              mc,nlie(mlev),llev(mmic,mlev)
c
c          nmic   = number of microphones [0,mmic]
c          posi   = spatial position (x/y/z) of each microphone [cm]
c          data   = data acquisition value from each microphone [Usec]
c          dist   = idem,                   converted in radius [cm]
c          lepo   = level position
c          nlie   = number of microphones per level
c          llev   = list of microphone numbers in each level
c
      common/spkdon/sbias(2),speed(2)
c
c          sbias  = time bias [Usec]
c          speed  = speed of sound in helium [cm/Usec]
c
      common/spkstr/string
      character*60 string
c
c          string = identification string for this data
c
+keep,spkevt.
c  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c
c  spark chamber event
c
      parameter(mpts=50,mtra=5)
      common/spkevt/numb,npts,coor(3,mpts),komi(2,mpts),mico(mpts),
     +                   ntra,ltra(2,mtra)
c
c          numb = event number
c
c          npts = number of points defined in space [0,mpts]
c          coor = spatial (x/y/z) coordinates of points [cm]
c          komi = from which 2 microphones/circles are coord. from
c          mico = code for microphone points [0=data,1=mc,2=average]
c
c          ntra = number of trajectories to be drawn [0,mtra]
c          ltra = list of pairs of points [1,npts] as trajectories
c
+patch,spkcod.   The SPK routines
+deck,spkprg.
c
c ======================================================================
      program spk
c ----------------------------------------------------------------------
c
c  purpose:  main program for the spark chamber display< etc..
c
c  history:  29.07.98  FC   created
c
c  comment:  preliminary..
c
c ----------------------------------------------------------------------
c
c  global declarations
c
+cde,spkios.
+cde,spkgeo.
+cde,spkdat.
+cde,spkevt.
c
c  local declarations
c
      real xy3(2,2),ds(mlev),dsmin(mlev)
      integer ii(mlev),jj(mlev),kk(mlev)
c
c  general initilizations
c
      call datime(ldate,ltime)
c
      lvers  =  20100
      lverd  = 020801
      lvert  =   1135
c
      lunin  = 5
      lunout = 6
      lunhpl = 7
      messag = 0
c
c  geometry
c
      platex =    30.00
      platey =    15.00
      spacez =     1.00
      thickz =     0.05
c
      nplate =    25
c
      bounds(1) = 45.00
      bounds(2) = 20.00
      bounds(3) = 20.00
c
c  time bias and speed of sound
c
      read(1,*) sbias(1),speed(1)
      read(1,*) sbias(2),speed(2)
c     write(8,*) '1',sbias,speed
c
c  microphone positions
c
      nmic = 0
      nlie(1) = 0
      nlie(2) = 0
c
1     if(nmic+1.gt.mmic) stop 1
      read(2,100,end=2) (posi(i,nmic+1),i=1,3)
100   format(3f10.0)
      nmic = nmic + 1
c     write(8,*) '2',nmic,(posi(i,nmic),i=1,3)
c
c     assume just 2 possible levels up to now
      if(posi(3,nmic).gt.0.0) then
         lep = 1
      else
         lep = 2
      endif
c
c     identify microphones in each level
      lepo(nmic) = lep
      nlie(lep) = nlie(lep) + 1
      llev(nlie(lep),lep) = nmic
c
c     get next microphone
      goto 1
c
c     end of microphones
2     if(nmic.lt.2) stop 2
c
c  plot initilizations
c
      open(unit=lunhpl,form='formatted',status='unknown')
      call spkini(-1,0)
c
c  identification string
c
      write(6,*) 'identification string?'
      read(5,101) string
101   format(a)
c
c  limit number of events to analyze
c
      write(6,*) 'number of events?'
      read(5,100) c1
      mev = abs(c1)
      if(mev.lt.1) mev = 1
      nev = 0
c
c  data (0) or monte-carlo (1)
c
      read(3,*) c3
      mc = c3
c
c  new event
c
3     nev = nev + 1
      if(nev.gt.mev) goto 4
c
c  microphone data
c
c     experimental
      if(mc.ne.1) then
         read(3,*,end=4) c1,c2,(data(i),i=1,nmic)
         npts = 0
         kpts = 0
         ntra = 0
c
c     monte-carlo, take known 'true' data as special points
      else
         read(3,*,end=4) c1,c2,(data(i),i=1,nmic),
     +                   coor(1,1),coor(1,2),coor(2,1),coor(2,2)
         coor(3,1) = posi(3,1)
         coor(3,2) = posi(3,2)
c
c        points
         npts = 2
         kpts = 2
         komi(1,1) = 0
         komi(2,1) = 0
         komi(1,2) = 0
         komi(2,2) = 0
         mico(1) = 1
         mico(2) = 1
c
c        trajectories
         ntra = 1
         ltra(1,1) = 1
         ltra(2,1) = 2
      endif
c
c     event info:  number & time
      numb = c1
      seco = c2 - 934944400.
c     write(8,*) '3',numb,seco
c     write(8,*) '4',(data(i),i=1,nmic)
c     if(mc.eq.1)
c    +write(8,*) '5',coor(1,1),coor(1,2),coor(2,1),coor(2,2)
c
c  convert to distance
c
      do i=1,nmic
         if(posi(3,i).gt.0.0) then
            j = 1
         else
            j = 2
         endif
         dist(i) = ( data(i)+sbias(j) ) * speed(j)
      enddo
c     write(8,*) '6',(dist(i),i=1,nmic)
c
c  draw available data before reconstruction
c
c     call spkplo(1.0)
c     call spkplo(0.2)
c     write(8,*) '61',(dist(i),i=1,nmic)
c
c  reconstruct event through the circle method
c
c     take only microphones from the same height
      do i=1,nmic-1
         do j=i+1,nmic
            if(posi(3,i).eq.posi(3,j)) then
c
c              find the intersection(s) of 2 circles
               call spkcir(posi(1,i),dist(i),posi(1,j),dist(j),
     +                     np,xy3)
c
c                 there is always at least one point, event if np=0
                  if(abs(xy3(1,1)).le.bounds(1).and.
     +               abs(xy3(2,1)).le.bounds(2)) then
                     npts = npts + 1
                     coor(1,npts) = xy3(1,1)
                     coor(2,npts) = xy3(2,1)
                     coor(3,npts) = posi(3,i)
                     komi(1,npts) = i
                     komi(2,npts) = j
                     mico(  npts) = 0
c                    write(8,*)
c    +                  '7',coor(1,npts),coor(2,npts),coor(3,npts)
c                    write(8,*)
c    +                  '8',npts,komi(1,npts),komi(2,npts)
                  endif
c
c              normal case of 2 points with ambiguity on which
               if(np.gt.1) then
                  if(abs(xy3(1,2)).le.bounds(1).and.
     +               abs(xy3(2,2)).le.bounds(2)) then
                     npts = npts + 1
                     coor(1,npts) = xy3(1,2)
                     coor(2,npts) = xy3(2,2)
                     coor(3,npts) = posi(3,i)
                     komi(1,npts) = i
                     komi(2,npts) = j
                     mico(  npts) = 0
c                    write(8,*)
c    +                  '9',coor(1,npts),coor(2,npts),coor(3,npts)
c                    write(8,*)
c    +                  '0',npts,komi(1,npts),komi(2,npts)
                  endif
               endif
c
            endif
         enddo
      enddo
c
c  try to find average intersection and resolve ambuiguities
c
c     local initializations
      do i=1,mlev
         dsmin(i) = 9999999.
         ds(i) = 0.0
         ii(i) = 0
         jj(i) = 0
         kk(i) = 0
      enddo
c
c     kpts was last MC 'true' point if any
c     lpts is the last data point, don't consider new averaged points
      lpts = npts
c
c     proceed level by level
      do il=1,mlev
      if(nlie(il).ge.3) then
c
c        assume no more than 3 micros and do the loops
cno      do i=1,npts-2
         do i=kpts+1,lpts-2
            i1 = komi(1,i)
            i2 = komi(2,i)
         do j=i,lpts-1
            j1 = komi(1,j)
            j2 = komi(2,j)
         do k=j,lpts
            k1 = komi(1,k)
            k2 = komi(2,k)
c
c        tried the following alternative, but something's wrong
c        do i=1,nlie(il)-2
c           i1 = komi(1,llev(i,il))
c           i2 = komi(2,llev(i,il))
c        do j=i,nlie(il)-1
c           j1 = komi(1,llev(j,il))
c           j2 = komi(2,llev(j,il))
c        do k=j,nlie(il)
c           k1 = komi(1,llev(k,il))
c           k2 = komi(2,llev(k,il))
c
c           write(8,*) 'a',i1,i2,j1,j2,k1,k2
c
c           this test becomes trivial if nlie-logic works
            if(posi(3,i1).eq.posi(3,j1).and.
     +         posi(3,i1).eq.posi(3,k1)) then
               lep = lepo(i1)
c
c           consider only different pair combinations of circles
            if((i1.ne.j1.or.i2.ne.j2).and.
     +         (i1.ne.k1.or.i2.ne.k2).and.
     +         (j1.ne.k1.or.j2.ne.k2)) then
c
c              compute some sum of distances between all
               ds(lep) =
     +            sqrt(
     +               (coor(1,i)-coor(1,j))**2+(coor(2,i)-coor(2,j))**2)+
     +            sqrt(
     +               (coor(1,i)-coor(1,k))**2+(coor(2,i)-coor(2,k))**2)+
     +            sqrt(
     +               (coor(1,j)-coor(1,k))**2+(coor(2,j)-coor(2,k))**2)
c
c              write(8,*) 'b',lep,dsmin(lep),ds(lep)
c              write(8,*) 'c',i,j,k
c
c              keep the combination with the lowest distance sum
               if(ds(lep).lt.dsmin(lep)) then
                  dsmin(lep) = ds(lep)
                  ii(lep) = i
                  jj(lep) = j
                  kk(lep) = k
               endif
c              write(8,*) 'd',lep,dsmin(lep),ds(lep)
            endif
            endif
c
c        end of combination
         enddo
         enddo
         enddo
c
c     end of combinations
      endif
      enddo
c
c     store averages, if any were found
      kpts = npts
      do i=1,mlev
         if(ii(i).gt.0) then
            npts = npts + 1
            coor(1,npts) =
     +         (coor(1,ii(i))+coor(1,jj(i))+coor(1,kk(i))) / 3.0
            coor(2,npts) =
     +         (coor(2,ii(i))+coor(2,jj(i))+coor(2,kk(i))) / 3.0
            coor(3,npts) = coor(3,ii(i))
            komi(1,npts) = -1
            komi(2,npts) = -1
            mico(  npts) = 2
         endif
      enddo
c
c     might be enough (assumes only 2 in fact) points for trajectory
      if(npts.ge.kpts+2) then
         ntra = ntra + 1
         ltra(1,ntra) = kpts + 1
         ltra(2,ntra) = kpts + 2
      endif
c
c  draw results after reconstruction
c
      if(npts.gt.0.or.ntra.gt.0) then
         call spkplo(1.0)
c        call spkplo(0.2)
      endif
c
c  next event
c
      goto 3
c
c  last event
c
4     continue
c
c  conclusions
c
      call spkini(+1,0)
      close(unit=lunhpl)
c
c  exit
c
      stop
      end
+deck,spkini.
c
c ======================================================================
      subroutine spkini(io,lterm)
c ----------------------------------------------------------------------
c
c  purpose:  initialize plotting for 'zut'
c
c  history:  29.07.98  FC   from zt1plo
c
c   inputs:  io     = -1/+1 for initialization/termination
c            lterm  = terminal type for 'hplot'
c
c  outputs:  none
c
c  comment:  meant as an example, since user may want larger /pawc/...
c
c ----------------------------------------------------------------------
c
c  global declarations
c
+cde,spkios.
c
c  local declarations
c
      parameter (mlim=20000)
      common/pawc/h(mlim)
c
c  initialize proper
c
      if(io.lt.0) then
         call hlimit(mlim)
         call houtpu(lunout)
         call hplint(lterm)
         if(lterm.ne.0) then
            call hplcap(+lunhpl)
         else
            call hplcap(-lunhpl)
         endif
         call igsa(0)
c
c  page size
c
         call hplsiz(21.0,27.8,' ')
c
c  graphics options
c
         call hplset('CSIZ',0.02)
         call hplset('VSIZ',0.01)
         call hplset('ASIZ',0.01)
         call hplset('TSIZ',0.01)
         call hplset('GSIZ',0.01)
         call hplset('XMGL',0.10)
         call hplset('XMGR',0.10)
         call hplset('YMGL',0.10)
         call hplset('YMGU',1.00)
         call hplset('XWIN',0.10)
         call hplset('YWIN',0.10)
         call hplset('XTIC',0.10)
         call hplset('YTIC',0.10)
         call hplset('YHTI',0.10)
         call hplset('YGTI',0.10)
         call hplset('XLAB',0.10)
         call hplset('YLAB',0.10)
         call hplset('XVAL',0.10)
c
c  terminate properly
c
      else
         call hplend
      endif
c
c  exit
c
      return
      end
+deck,spklin.
c
c ======================================================================
      subroutine spklin(xx,yy)
c ----------------------------------------------------------------------
c
c  purpose:  plot one line
c
c  history:  29.07.98  FC   from zt3lin
c
c   inputs:  xx(2)  = x-co-ordinates
c            yy(2)  = y-co-ordinates
c
c  outputs:  none
c
c ----------------------------------------------------------------------
c
c  local declarations
c
      real x(2),y(2),xx(2),yy(2)
c
c  check sequence
c
      if(xx(1).gt.xx(2)) then
         x(1) = xx(2)
         x(2) = xx(1)
         y(1) = yy(2)
         y(2) = yy(1)
      else
         x(1) = xx(1)
         x(2) = xx(2)
         y(1) = yy(1)
         y(2) = yy(2)
      endif
c
c  make line
c
      call hpline(x,y,2,' ')
c
c  exit
c
      return
      end
+deck,spkplo.
c
c ======================================================================
      subroutine spkplo(zoomfa)
c ----------------------------------------------------------------------
c
c  purpose:  plot event graphically
c
c  history:  29.07.98  FC   from zt3plo
c
c   inputs:  zoomfa = zooming factor, ]1, [ = larger frame in fact
c                                     ]0,1[ = real zooming
c                                     =0=1  = no zooming
c                                     <0    = same but energy vectors
c
c  outputs:  none
c
c ----------------------------------------------------------------------
c
c  global declarations
c
+cde,spkgeo.
c
c  histogram definitions
c
      zoom  = abs(zoomfa)
      if(zoom.eq.0.0) zoom = 1.0
c
      call hbook2(1,'yz', 2,-bounds(2)*zoom,+bounds(2)*zoom,
     +                    2,-bounds(3)*zoom,+bounds(3)*zoom, 0.0)
      call hbook2(2,'xz', 2,-bounds(1)*zoom,+bounds(1)*zoom,
     +                    2,-bounds(3)*zoom,+bounds(3)*zoom, 0.0)
      call hbook2(3,'xy', 2,-bounds(1)*zoom,+bounds(1)*zoom,
     +                    2,-bounds(2)*zoom,+bounds(2)*zoom, 0.0)
c
c  zeroth zone
c
c     write(8,*) '6a',zoomfa
      call hplzon(2,3,1,' ')
      call spkzon(1,0,0,zoomfa)
c
c  first zone
c
      call hplzon(2,3,2,'S')
      call spkzon(1,2,3,zoomfa)
      call hdelet(1)
c
c  second zone
c
      call hplzon(1,3,2,'S')
      call spkzon(2,1,3,zoomfa)
      call hdelet(2)
c
c  third zone
c
      call hplzon(1,3,3,'S')
      call spkzon(3,1,2,zoomfa)
      call hdelet(3)
c
c  exit
c
      return
      end
+deck,spkzon.
c
c ======================================================================
      subroutine spkzon(no,ix,iy,zoomfa)
c ----------------------------------------------------------------------
c
c  purpose:  plot one zone of event
c
c  history:  29.07.98  FC   from zt3zon
c
c   inputs:  no     = histogram/zone number
c            ix     = coordinate for horizontal axis
c            iy     = .          .   vertical   .
c            zoomfa = zooming factor, ]1, [ = larger frame in fact
c                                     ]0,1[ = real zooming
c                                     <0    = same but energy vectors
c
c  outputs:  none
c
c ----------------------------------------------------------------------
c
c  global declarations
c
+cde,spkios.
+cde,spkgeo.
+cde,spkdat.
+cde,spkevt.
c
c  local declarations
c
      real xx(2),yy(2)
      character*80 line
c
c  cut parameters for the projections
c
      real xcut,ycut,zcut
      save xcut,ycut,zcut
      data xcut,ycut,zcut/0.0,0.0,0.0/
c
c  zoom factor
c
      zoom  = abs(zoomfa)
      if(zoom.eq.0.0) zoom = 1.0
c
c  plot frame
c
      call hplot(no,' ','hist',0)
c
c  patterns
c
c     call isfais(3)
      call igset('FAIS',0.0)
      call igset('BORD',1.0)
c
c  common to all
c
      z0 = - float(nplate-1) * spacez / 2.0
      dz = thickz / 2.0
c
c  description
c
      if(ix*iy.eq.0) then
c
c        code
         write(line,202) ldate,ltime,lvers,lverd,lvert
202      format('date=',i8,'  time=',i4,'  vers=',i5,' (',i6,i5,')')
         call igtext((-0.96*bounds(2))*zoom,
     +               (+0.92*bounds(3))*zoom,line,
     +               (+0.06*bounds(2))*0.5*zoom,0.0,'L')
c
c        symbols
         call hplsym((-0.94*bounds(2))*zoom,
     +               (+0.82*bounds(3))*zoom,
     +               1,31,0.2,' ')
         call igtext((-0.90*bounds(2))*zoom,
     +               (+0.80*bounds(3))*zoom,'origin',
     +               (+0.06*bounds(2))*0.5*zoom,0.0,'L')
c
         call hplsym((-0.94*bounds(2))*zoom,
     +               (+0.76*bounds(3))*zoom,
     +               1,24,0.2,' ')
         call igtext((-0.90*bounds(2))*zoom,
     +               (+0.74*bounds(3))*zoom,'microphones z"g#0',
     +               (+0.06*bounds(2))*0.5*zoom,0.0,'L')
c
         call hplsym((-0.94*bounds(2))*zoom,
     +               (+0.70*bounds(3))*zoom,
     +               1,25,0.2,' ')
         call igtext((-0.90*bounds(2))*zoom,
     +               (+0.68*bounds(3))*zoom,'microphones z"l#0',
     +               (+0.06*bounds(2))*0.5*zoom,0.0,'L')
c
         call hplsym((-0.94*bounds(2))*zoom,
     +               (+0.64*bounds(3))*zoom,
     +               1,20,0.2,' ')
         call igtext((-0.90*bounds(2))*zoom,
     +               (+0.62*bounds(3))*zoom,'hits z"g#0',
     +               (+0.06*bounds(2))*0.5*zoom,0.0,'L')
c
         call hplsym((-0.94*bounds(2))*zoom,
     +               (+0.58*bounds(3))*zoom,
     +               1,21,0.2,' ')
         call igtext((-0.90*bounds(2))*zoom,
     +               (+0.56*bounds(3))*zoom,'hits z"l#0',
     +               (+0.06*bounds(2))*0.5*zoom,0.0,'L')
c
         call hplsym((-0.94*bounds(2))*zoom,
     +               (+0.52*bounds(3))*zoom,
     +               1,22,0.2,' ')
         call igtext((-0.90*bounds(2))*zoom,
     +               (+0.50*bounds(3))*zoom,'average z"g#0',
     +               (+0.06*bounds(2))*0.5*zoom,0.0,'L')
c
         call hplsym((-0.94*bounds(2))*zoom,
     +               (+0.46*bounds(3))*zoom,
     +               1,23,0.2,' ')
         call igtext((-0.90*bounds(2))*zoom,
     +               (+0.44*bounds(3))*zoom,'average z"l#0',
     +               (+0.06*bounds(2))*0.5*zoom,0.0,'L')
c
         call hplsym((-0.94*bounds(2))*zoom,
     +               (+0.40*bounds(3))*zoom,
     +               1,28,0.2,' ')
         call igtext((-0.90*bounds(2))*zoom,
     +               (+0.38*bounds(3))*zoom,'Monte-Carlo',
     +               (+0.06*bounds(2))*0.5*zoom,0.0,'L')
c
         call igtext((-0.90*bounds(2))*zoom,
     +               (+0.32*bounds(3))*zoom,'full"j# z"g#0',
     +               (+0.06*bounds(2))*0.5*zoom,0.0,'L')
c
         call igtext((-0.90*bounds(2))*zoom,
     +               (+0.26*bounds(3))*zoom,'dashed"j# z"l#0',
     +               (+0.06*bounds(2))*0.5*zoom,0.0,'L')
c
         call igtext((-0.90*bounds(2))*zoom,
     +               (+0.20*bounds(3))*zoom,'dotted"j# MC true',
     +               (+0.06*bounds(2))*0.5*zoom,0.0,'L')
c
c        data & event & info
         write(line,206) string
206      format('ID =',1x,a60)
         call igtext((-0.96*bounds(2))*zoom,
     +               (+0.10*bounds(3))*zoom,line,
     +               (+0.06*bounds(2))*0.5*zoom,0.0,'L')
         write(line,203) mc,numb
203      format('MC,event=',i2,i6)
         call igtext((-0.96*bounds(2))*zoom,
     +               (+0.04*bounds(3))*zoom,line,
     +               (+0.06*bounds(2))*0.5*zoom,0.0,'L')
c
         if(mc.ne.1) then
         write(line,204) sbias
204      format('sbias=',2f8.3)
         call igtext((-0.96*bounds(2))*zoom,
     +               (-0.06*bounds(3))*zoom,line,
     +               (+0.06*bounds(2))*0.5*zoom,0.0,'L')
         write(line,205) speed
205      format('speed=',2f8.3)
         call igtext((-0.96*bounds(2))*zoom,
     +               (-0.12*bounds(3))*zoom,line,
     +               (+0.06*bounds(2))*0.5*zoom,0.0,'L')
         endif
         return
c
c  plot detector (yz,xz,xy)
c
c  yz-view
c
      else if(no.eq.1) then
         call igtext((-0.96*bounds(2))*zoom,
     +               (+0.88*bounds(3))*zoom,'yz',
     +               (+0.06*bounds(2))*zoom,0.0,'L')
c
         do ipl=1,nplate
            x = platey
            y = z0 + float(ipl-1) * spacez
            call igbox(-x,+x,y-dz,y+dz)
         enddo
c
c  xz-view
c
      else if(no.eq.2) then
         call igtext((-0.96*bounds(1))*zoom,
     +               (+0.88*bounds(3))*zoom,'xz',
     +               (+0.06*bounds(2))*zoom,0.0,'L')
c
         do ipl=1,nplate
            x = platex
            y = z0 + float(ipl-1) * spacez
            call igbox(-x,+x,y-dz,y+dz)
         enddo
c
c  xy-view
c
      else if(no.eq.3) then
         call igtext((-0.96*bounds(1))*zoom,
     +               (+0.88*bounds(2))*zoom,'xy',
     +               (+0.06*bounds(2))*zoom,0.0,'L')
c
         do ipl=1,nplate
            x = platex
            y = platey
            call igbox(-x,+x,-y,+y)
         enddo
c
      endif
c
c  start of data
c
      if(zoomfa.ge.0.0) then
c
c  origin
c
         call hplsym(0.0,0.0,1,31,0.2,' ')
c
c  plot microphone information
c
         if(nmic.gt.0) then
            do i=1,nmic
               if(posi(3,i).gt.0.0) then
                  isym = 24
               else
                  isym = 25
               endif
               call hplsym(posi(ix,i),posi(iy,i),1,isym,0.2,' ')
            enddo
         endif
c
c  plot circles around microphones
c
         if(ix.ne.3.and.iy.ne.3) then
            do i=1,nmic
               if(posi(3,i).gt.0.0) then
                  call igset('LTYP',1.0)
               else
                  call igset('LTYP',2.0)
               endif
               call igarc(posi(ix,i),posi(iy,i),dist(i),dist(i),0.,0.)
            enddo
         endif
c
c  plot point information
c
         if(npts.gt.0) then
            do i=1,npts
c
c              true MC points
	       if(mico(i).eq.1) then
		  if(coor(3,i).gt.0.0) then
		     isym = 28
		  else
		     isym = 28
		  endif
c
c              average points
	       else if(mico(i).eq.2) then
		  if(coor(3,i).gt.0.0) then
		     isym = 22
		  else
		     isym = 23
		  endif
c
c              normal data points
	       else
		  if(coor(3,i).gt.0.0) then
		     isym = 20
		  else
		     isym = 21
		  endif
	       endif
c
c              draw the point
	       call hplsym(coor(ix,i),coor(iy,i),1,isym,0.2,' ')
	    enddo
c
c  plot trajectory information
c
	    if(ntra.gt.0) then
	       do i=1,ntra
c
c                 the 2 known points of this trajectory
		  n1 = ltra(1,i)
		  n2 = ltra(2,i)
c                 write(8,*) 'e',i,n1,n2
c
c                 coordinates of the valid points
		  if(n1.ge.1.and.n1.le.npts.and.
     +               n2.ge.1.and.n2.le.npts.and.
     +               n1.ne.n2) then
		     x1 = coor(ix,n1)
		     y1 = coor(iy,n1)
		     x2 = coor(ix,n2)
                     y2 = coor(iy,n2)
c
c                    extrapolate to drawing boundaries
                     if(x1.ne.x2.or.y1.ne.y2) then
c
c                       depends which way around
                        if(x1.ne.x2) then
                           if(x1.lt.x2) then
                              a = (y2-y1)/(x2-x1)
                           else
                              a = (y1-y2)/(x1-x2)
                           endif
                           b = y1 - a * x1
                           xx(1) = -bounds(ix)
                           yy(1) = a * xx(1) + b
                           xx(2) = +bounds(ix)
                           yy(2) = a * xx(2) + b
                        else
                           xx(1) = x1
                           xx(2) = x2
                           yy(1) = -bounds(iy)
                           yy(2) = +bounds(iy)
                        endif
c
c                       do jkl=1,npts
c                       write(8,*) 'f',jkl,mico(jkl),(coor(m,jkl),m=1,3)
c                       enddo
c
c                       line type of the trajectory
                        if(mico(n1).ne.1) then
                           call hplset('DMOD',1.0)
c                          write(8,*) 'g',n1,mico(n1)
                        else
                           call hplset('DMOD',4.0)
c                          write(8,*) 'h',n1,mico(n1)
                        endif
c
c                       draw the trajectory
                        call spklin(xx,yy)
                     endif
                  endif
               enddo
c
c  end of trajectories
c
            endif
c
c  end of points
c
         endif
c
c  end of data
c
      endif
c
c  exit
c
      return
      end
+deck,spkcir.
c
c ======================================================================
      subroutine spkcir(xy1,r1,xy2,r2,np,xy3)
c ----------------------------------------------------------------------
c
c  purpose:  find instersections of 2 circles
c
c  history:  23.08.99  fc   from scratch
c
c   inputs:  xy1(2)    center of 1st circle
c            r1        radius .  .   .
c            xy2(2)    center .  2nd .
c            r2        radius .  .   .
c
c  outputs:  np        number of intersection points [0,1,2][-1:error]
c            xy3(2,2)  coord. .  .            .      (2,np)
c
c ----------------------------------------------------------------------
c
c  declarations
c
      real xy1(2),xy2(2),xy3(2,2)
c
c  initializations
c
      np = 0
      do i=1,2
         do j=1,2
            xy3(j,i) = 0.0
         enddo
      enddo
c
c  distance between both centers
c
      d = sqrt( (xy1(1)-xy2(1))**2 + (xy1(2)-xy2(2))**2 )
c     write(8,*) 'i',r1,r2,d
c
      if(d.le.0.0) then
         np = -1
         return
      endif
c
c  changes of coordinates
c
c     rotate around 1st point to place 2nd on x-axis
c     only the distance d is relevant since:
c     (x1,y1) & (x2,y2) --> (0,0) & (d,0) in new system
c
c     rotate back
      cost = (xy2(1)-xy1(1)) / d
      sint = (xy2(2)-xy1(2)) / d
c
c  case 1:  circles don't touch
c
      if(d.gt.r1+r2) then
         np = 0
c        take something anyway:  halfway, closest approach
         xp = ( r1 + (d-r2) ) / 2.0
         yp = 0.0
c        write(8,*) 'j',xp,yp,np
c
c  case 2:  circles just touch
c
      else if(d.eq.r1+r2) then
         np = 1
c        take this point
         xp = r1
         yp = 0.0
c        write(8,*) 'k',xp,yp,np
c
c  case 3:  regular intersection, there will be 2 solutions
c
      else
         np = 2
c        solution of equation, +-yp are the two cases
         xp = ( r1**2 - ( r2**2 - d**2 ) ) / ( 2.0 * d )
         yp = sqrt( r1**2 - xp**2 )
c        write(8,*) 'l',xp,yp,np
      endif
c
c  convert back
c
c     always at least one case is considered
         xy3(1,1) = xp * cost - yp * sint + xy1(1)
         xy3(2,1) = xp * sint + yp * cost + xy1(2)
c        write(8,*) 'm',xy3(1,1),xy3(2,1)
c
c     the normal case has 2 points
      if(np.gt.1) then
         xy3(1,2) = xp * cost + yp * sint + xy1(1)
         xy3(2,2) = xp * sint - yp * cost + xy1(2)
c        write(8,*) 'n',xy3(1,2),xy3(2,2)
      endif
c
c  exit
c
      return
      end