1 SUBROUTINE PXCONE(MODE,NTRAK,ITKDM,PTRAK,CONER,EPSLON,OVLIM,
2 + MXJET,NJET,PJET,IPASS,IJMUL,IERR)
3 *.*********************************************************
8 *.********** Pre Release Version 26.2.93
10 *. Driver for the Cone Jet finding algorithm of L.A. del Pozo.
11 *. Based on algorithm from D.E. Soper.
12 *. Finds jets inside cone of half angle CONER with energy > EPSLON.
13 *. Jets which receive more than a fraction OVLIM of their energy from
14 *. overlaps with other jets are excluded.
15 *. Output jets are ordered in energy.
16 *. If MODE.EQ.2 momenta are stored as (eta,phi,<empty>,pt)
19 *. INTEGER ITKDM,MXTRK
20 *. PARAMETER (ITKDM=4.or.more,MXTRK=1.or.more)
21 *. INTEGER MXJET, MXTRAK, MXPROT
22 *. PARAMETER (MXJET=10,MXTRAK=500,MXPROT=500)
23 *. INTEGER IPASS (MXTRAK),IJMUL (MXJET)
24 *. INTEGER NTRAK,NJET,IERR,MODE
25 *. DOUBLE PRECISION PTRAK (ITKDM,MXTRK),PJET (5,MXJET)
26 *. DOUBLE PRECISION CONER, EPSLON, OVLIM
27 *. NTRAK = 1.to.MXTRAK
31 *. CALL PXCONE (MODE,NTRAK,ITKDM,PTRAK,CONER,EPSLON,OVLIM,MXJET,
32 *. + NJET,PJET,IPASS,IJMUL,IERR)
34 *. INPUT : MODE 1=>e+e-, 2=>hadron-hadron
35 *. INPUT : NTRAK Number of particles
36 *. INPUT : ITKDM First dimension of PTRAK array
37 *. INPUT : PTRAK Array of particle 4-momenta (Px,Py,Pz,E)
38 *. INPUT : CONER Cone size (half angle) in radians
39 *. INPUT : EPSLON Minimum Jet energy (GeV)
40 *. INPUT : OVLIM Maximum fraction of overlap energy in a jet
41 *. INPUT : MXJET Maximum possible number of jets
42 *. OUTPUT : NJET Number of jets found
43 *. OUTPUT : PJET 5-vectors of jets
44 *. OUTPUT : IPASS(k) Particle k belongs to jet number IPASS(k)
45 *. IPASS = -1 if not assosciated to a jet
46 *. OUTPUT : IJMUL(i) Jet i contains IJMUL(i) particles
47 *. OUTPUT : IERR = 0 if all is OK ; = -1 otherwise
49 *. CALLS : PXSEAR, PXSAME, PXNEW, PXTRY, PXORD, PXUVEC, PXOLAP
52 *. AUTHOR : L.A. del Pozo
53 *. CREATED : 26-Feb-93
54 *. LAST MOD : 2-Mar-93
57 *. 2-Jan-97: M Wobisch - fix bug concerning COS2R in eta phi mode
58 *. 4-Apr-93: M H Seymour - Change 2d arrays to 1d in PXTRY & PXNEW
59 *. 2-Apr-93: M H Seymour - Major changes to add boost-invariant mode
60 *. 1-Apr-93: M H Seymour - Increase all array sizes
61 *. 30-Mar-93: M H Seymour - Change all REAL variables to DOUBLE PRECISION
62 *. 30-Mar-93: M H Seymour - Change OVLIM into an input parameter
63 *. 2-Mar-93: L A del Pozo - Fix Bugs in PXOLAP
64 *. 1-Mar-93: L A del Pozo - Remove Cern library routine calls
65 *. 1-Mar-93: L A del Pozo - Add Print out of welcome and R and Epsilon
67 *.*********************************************************
70 INTEGER ITKDM,MXJET,NTRAK,NJET,IERR,MODE
71 INTEGER IPASS (NTRAK),IJMUL (MXJET)
72 DOUBLE PRECISION PTRAK (ITKDM,NTRAK),PJET (5,MXJET),
73 + CONER, EPSLON, OVLIM
75 INTEGER MXPROT, MXTRAK
76 PARAMETER (MXPROT=5000, MXTRAK=5000)
77 DOUBLE PRECISION PP(4,MXTRAK), PU(3,MXTRAK), PJ(4,MXPROT)
78 LOGICAL JETLIS(MXPROT,MXTRAK)
79 *** Used in the routine.
80 DOUBLE PRECISION COSR,COS2R, VSEED(3), VEC1(3), VEC2(3),PTSQ,PPSQ,
86 INTEGER I,J,N,MU,N1,N2, ITERR
88 DOUBLE PRECISION ROLD, EPSOLD, OVOLD
89 SAVE NCALL,NPRINT,ROLD, EPSOLD, OVOLD
90 DATA NCALL,NPRINT /0,0/
91 DATA ROLD,EPSOLD,OVOLD/0.,0.,0./
94 c***************************************
96 c***************************************
108 *** Print welcome and Jetfinder parameters
109 IF((CONER.NE.ROLD .OR. EPSLON.NE.EPSOLD .OR. OVLIM.NE.OVOLD)
110 + .AND. NPRINT.LE.10) THEN
112 WRITE (6,*) ' *********** PXCONE: Cone Jet-finder ***********'
113 WRITE (6,*) ' Written by Luis Del Pozo of OPAL'
114 WRITE (6,*) ' Modified for eta-phi by Mike Seymour'
115 WRITE(6,1000)' Cone Size R = ',CONER,' Radians'
116 WRITE(6,1001)' Min Jet energy Epsilon = ',EPSLON,' GeV'
117 WRITE(6,1002)' Overlap fraction parameter = ',OVLIM
118 WRITE (6,*) ' ***********************************************'
120 IF (RSEP .lt. 1.999) THEN
122 WRITE (6,*) ' ******************************************'
123 WRITE (6,*) ' ******************************************'
124 WRITE(6,*) ' M Wobisch: private change !!!!!!!!!!!! '
125 WRITE(6,*) ' Rsep is set to ',RSEP
126 WRITE(6,*) ' this is ONLY meaningful in a NLO calculation'
127 WRITE(6,*) ' ------------------------ '
128 WRITE(6,*) ' please check what you''re doing!!'
129 WRITE(6,*) ' or ask: Markus.Wobisch@desy.de --'
130 WRITE (6,*) ' ******************************************'
131 WRITE (6,*) ' ******************************************'
132 WRITE (6,*) ' ******************************************'
139 1000 FORMAT(A18,F5.2,A10)
140 1001 FORMAT(A29,F5.2,A5)
141 1002 FORMAT(A33,F5.2)
148 *** Copy calling array PTRAK to internal array PP(4,NTRAK)
150 IF (NTRAK .GT. MXTRAK) THEN
151 WRITE (6,*) ' PXCONE: Ntrak too large: ',NTRAK
162 *** Converting to eta,phi,pt if necessary
164 PTSQ=PTRAK(1,I)**2+PTRAK(2,I)**2
165 PPSQ=(SQRT(PTSQ+PTRAK(3,I)**2)+ABS(PTRAK(3,I)))**2
166 IF (PTSQ.LE.4.25E-18*PPSQ) THEN
169 PP(1,I)=0.5*LOG(PPSQ/PTSQ)
171 PP(1,I)=SIGN(PP(1,I),PTRAK(3,I))
175 PP(2,I)=ATAN2(PTRAK(2,I),PTRAK(1,I))
185 *** Zero output variables
190 JETLIS(J,I) = .FALSE.
193 CALL PXZERV(4*MXPROT,PJ)
194 CALL PXZERI(MXJET,IJMUL)
200 *** Purely for convenience, work in terms of 1-R**2
202 cMW -- select Rsep: 1-(Rsep*CONER)**2
203 COS2R = 1-(RSEP*CONER)**2
204 cORIGINAL COS2R = 1-(2*CONER)**2
208 CALL PXUVEC(NTRAK,PP,PU,IERR)
209 IF (IERR .NE. 0) RETURN
211 *** Look for jets using particle diretions as seed axes
217 CALL PXSEAR(MODE,COSR,NTRAK,PU,PP,VSEED,
218 & NJET,JETLIS,PJ,UNSTBL,IERR)
219 IF (IERR .NE. 0) RETURN
222 cMW - for Rsep=1 goto 145
225 *** Now look between all pairs of jets as seed axes.
230 IF (MODE.NE.2) CALL PXNORV(3,VEC1,VEC1,ITERR)
231 DO 150 N2 = N1+1,NJET
235 IF (MODE.NE.2) CALL PXNORV(3,VEC2,VEC2,ITERR)
236 CALL PXADDV(3,VEC1,VEC2,VSEED,ITERR)
238 CALL PXNORV(3,VSEED,VSEED,ITERR)
243 C---ONLY BOTHER IF THEY ARE BETWEEN 1 AND 2 CONE RADII APART
245 COSVAL=VEC1(1)*VEC2(1)+VEC1(2)*VEC2(2)+VEC1(3)*VEC2(3)
247 IF (ABS(VEC1(1)).GE.20.OR.ABS(VEC2(1)).GE.20) THEN
251 + ((VEC1(1)-VEC2(1))**2+PXMDPI(VEC1(2)-VEC2(2))**2)
255 IF (COSVAL.LE.COSR.AND.COSVAL.GE.COS2R)
256 + CALL PXSEAR(MODE,COSR,NTRAK,PU,PP,VSEED,NJET,
257 + JETLIS,PJ,UNSTBL,IERR)
258 c CALL PXSEAR(MODE,COSR,NTRAK,PU,PP,VSEED,NJET,
259 c + JETLIS,PJ,UNSTBL,IERR)
260 IF (IERR .NE. 0) RETURN
265 WRITE (6,*) ' PXCONE: Too many iterations to find a proto-jet'
270 *** Now put the jet list into order by jet energy, eliminating jets
271 *** with energy less than EPSLON.
272 CALL PXORD(EPSLON,NJET,NTRAK,JETLIS,PJ)
274 *** Take care of jet overlaps
275 CALL PXOLAP(MODE,NJET,NTRAK,JETLIS,PJ,PP,OVLIM)
277 *** Order jets again as some have been eliminated, or lost energy.
278 CALL PXORD(EPSLON,NJET,NTRAK,JETLIS,PJ)
280 *** All done!, Copy output into output arrays
281 IF (NJET .GT. MXJET) THEN
282 WRITE (6,*) ' PXCONE: Found more than MXJET jets'
294 PJET(1,I)=PJ(4,I)*COS(PJ(2,I))
295 PJET(2,I)=PJ(4,I)*SIN(PJ(2,I))
296 PJET(3,I)=PJ(4,I)*SINH(PJ(1,I))
297 PJET(4,I)=PJ(4,I)*COSH(PJ(1,I))
303 IF (JETLIS(J,I)) THEN
311 *CMZ : 1.06/00 28/02/94 15.44.44 by P. Schleper
313 C-----------------------------------------------------------------------
314 SUBROUTINE PXNORV(N,A,B,ITERR)
316 DOUBLE PRECISION A(N),B(N),C
327 *CMZ : 2.00/00 10/01/95 10.17.57 by P. Schleper
328 *CMZ : 1.06/00 15/03/94 12.17.46 by P. Schleper
332 SUBROUTINE PXOLAP(MODE,NJET,NTRAK,JETLIS,PJ,PP,OVLIM)
334 *** Looks for particles assigned to more than 1 jet, and reassigns them
335 *** If more than a fraction OVLIM of a jet's energy is contained in
336 *** higher energy jets, that jet is neglected.
337 *** Particles assigned to the jet closest in angle (a la CDF, Snowmass).
339 INTEGER MXTRAK, MXPROT
340 PARAMETER (MXTRAK=5000,MXPROT=5000)
341 INTEGER NJET, NTRAK, MODE
342 LOGICAL JETLIS(MXPROT,MXTRAK)
343 DOUBLE PRECISION PJ(4,MXPROT),PP(4,MXTRAK),PXMDPI
346 DOUBLE PRECISION EOVER
347 DOUBLE PRECISION OVLIM
348 INTEGER ITERR, IJMIN, IJET(MXPROT), NJ
349 DOUBLE PRECISION VEC1(3), VEC2(3), COST, THET, THMIN
352 IF (NJET.LE.1) RETURN
353 *** Look for jets with large overlaps with higher energy jets.
355 *** Find overlap energy between jets I and all higher energy jets.
360 IF (JETLIS(I,N).AND.JETLIS(J,N)) THEN
365 EOVER = EOVER + PP(4,N)
368 *** Is the fraction of energy shared larger than OVLIM?
369 IF (EOVER.GT.OVLIM*PJ(4,I)) THEN
370 *** De-assign all particles from Jet I
372 JETLIS(I,N) = .FALSE.
376 *** Now there are no big overlaps, assign every particle in
377 *** more than 1 jet to the closet jet.
378 *** Any particles now in more than 1 jet are assigned to the CLOSET
389 *** Particle in > 1 jet - calc angles...
395 VEC2(1)=PJ(1,IJET(J))
396 VEC2(2)=PJ(2,IJET(J))
397 VEC2(3)=PJ(3,IJET(J))
399 CALL PXANG3(VEC1,VEC2,COST,THET,ITERR)
401 THET=(VEC1(1)-VEC2(1))**2+PXMDPI(VEC1(2)-VEC2(2))**2
406 ELSEIF (THET .LT. THMIN) THEN
411 *** Assign track to IJMIN
413 JETLIS(J,I) = .FALSE.
415 JETLIS(IJMIN,I)=.TRUE.
424 IF( JETLIS(I,N) ) THEN
427 PJ(MU,I) = PJ(MU,I) + PP(MU,N)
431 + + PP(4,N)/(PP(4,N)+PJ(4,I))*(PP(1,N)-PJ(1,I))
433 + + PP(4,N)/(PP(4,N)+PJ(4,I))*PXMDPI(PP(2,N)-PJ(2,I))
434 PJ(4,I)=PJ(4,I)+PP(4,N)
441 *CMZ : 2.00/00 10/01/95 10.17.57 by P. Schleper
442 *CMZ : 1.06/00 14/03/94 15.37.45 by P. Schleper
446 SUBROUTINE PXORD(EPSLON,NJET,NTRAK,JETLIS,PJ)
448 *** Routine to put jets into order and eliminate tose less than EPSLON
450 INTEGER MXTRAK,MXPROT
451 PARAMETER (MXTRAK=5000,MXPROT=5000)
452 INTEGER I, J, INDEX(MXPROT)
453 DOUBLE PRECISION PTEMP(4,MXPROT), ELIST(MXPROT)
455 LOGICAL JETLIS(MXPROT,MXTRAK)
456 LOGICAL LOGTMP(MXPROT,MXTRAK)
457 DOUBLE PRECISION EPSLON,PJ(4,MXPROT)
458 *** Puts jets in order of energy: 1 = highest energy etc.
459 *** Then Eliminate jets with energy below EPSLON
461 *** Copy input arrays.
467 LOGTMP(I,J)=JETLIS(I,J)
473 *** Sort the energies...
474 CALL PXSORV(NJET,ELIST,INDEX,'I')
475 *** Fill PJ and JETLIS according to sort ( sort is in ascending order!!)
478 PJ(J,I)=PTEMP(J,INDEX(NJET+1-I))
481 JETLIS(I,J)=LOGTMP(INDEX(NJET+1-I),J)
484 ** Jets are now in order
485 *** Now eliminate jets with less than Epsilon energy
487 IF (PJ(4,I) .LT. EPSLON) THEN
495 ********************************************************************
496 *CMZ : 2.00/00 10/01/95 10.17.57 by P. Schleper
497 *CMZ : 1.06/00 14/03/94 15.37.44 by P. Schleper
500 SUBROUTINE PXSEAR(MODE,COSR,NTRAK,PU,PP,VSEED,NJET,
501 + JETLIS,PJ,UNSTBL,IERR)
504 INTEGER MXTRAK, MXPROT
505 PARAMETER (MXTRAK=5000,MXPROT=5000)
506 INTEGER NTRAK, IERR, MODE
507 DOUBLE PRECISION COSR,PU(3,MXTRAK),PP(4,MXTRAK),VSEED(3)
509 LOGICAL JETLIS(MXPROT,MXTRAK)
511 DOUBLE PRECISION PJ(4,MXPROT)
512 *** Using VSEED as a trial axis , look for a stable jet.
513 *** Check stable jets against those already found and add to PJ.
514 *** Will try up to MXITER iterations to get a stable set of particles
517 LOGICAL PXSAME,PXNEW,OK
518 LOGICAL NEWLIS(MXTRAK),OLDLIS(MXTRAK)
519 DOUBLE PRECISION OAXIS(3),NAXIS(3),PNEW(4)
521 PARAMETER(MXITER = 30)
524 OAXIS(MU) = VSEED(MU)
529 DO 120 ITER = 1,MXITER
530 CALL PXTRY(MODE,COSR,NTRAK,PU,PP,OAXIS,NAXIS,PNEW,NEWLIS,OK)
531 *** Return immediately if there were no particles in the cone.
535 IF(PXSAME(NEWLIS,OLDLIS,NTRAK)) THEN
536 *** We have a stable jet.
537 IF (PXNEW(NEWLIS,JETLIS,NTRAK,NJET)) THEN
538 *** And the jet is a new one. So add it to our arrays.
539 *** Check arrays are big anough...
540 IF (NJET .EQ. MXPROT) THEN
541 WRITE (6,*) ' PXCONE: Found more than MXPROT proto-jets'
547 JETLIS(NJET,N) = NEWLIS(N)
555 *** The jet was not stable, so we iterate again
566 *CMZ : 1.06/00 28/02/94 15.44.44 by P. Schleper
568 C-----------------------------------------------------------------------
569 SUBROUTINE PXSORV(N,A,K,OPT)
570 C Sort A(N) into ascending order
571 C OPT = 'I' : return index array K only
572 C OTHERWISE : return sorted A and index array K
573 C-----------------------------------------------------------------------
575 PARAMETER (NMAX=5000)
577 * INTEGER N,I,J,K(N),IL(NMAX),IR(NMAX)
579 INTEGER N,I,J,K(NMAX),IL(NMAX),IR(NMAX)
582 * DOUBLE PRECISION A(N),B(NMAX)
583 DOUBLE PRECISION A(NMAX),B(NMAX)
585 IF (N.GT.NMAX) STOP 'Sorry, not enough room in Mike''s PXSORV'
592 2 IF(A(I).GT.A(J)) GO TO 5
593 3 IF(IL(J).EQ.0) GO TO 4
599 5 IF(IR(J).LE.0) GO TO 6
609 8 IF(IL(J).GT.0) GO TO 20
618 30 IF(OPT.EQ.'I') RETURN
623 *********************************************************************
624 *CMZ : 2.00/00 10/01/95 10.17.57 by P. Schleper
625 *CMZ : 1.06/00 14/03/94 15.37.44 by P. Schleper
629 SUBROUTINE PXTRY(MODE,COSR,NTRAK,PU,PP,OAXIS,NAXIS,
634 PARAMETER (MXTRAK=5000)
636 *** Note that although PU and PP are assumed to be 2d arrays, they
637 *** are used as 1d in this routine for efficiency
638 DOUBLE PRECISION COSR,PU(3*MXTRAK),PP(4*MXTRAK),OAXIS(3),PXMDPI
640 LOGICAL NEWLIS(MXTRAK)
641 DOUBLE PRECISION NAXIS(3),PNEW(4)
642 *** Finds all particles in cone of size COSR about OAXIS direction.
643 *** Calculates 4-momentum sum of all particles in cone (PNEW) , and
644 *** returns this as new jet axis NAXIS (Both unit Vectors)
646 DOUBLE PRECISION COSVAL,NORMSQ,NORM
660 COSVAL=COSVAL+OAXIS(MU)*PU(MU+NPU)
663 IF (ABS(PU(1+NPU)).GE.20.OR.ABS(OAXIS(1)).GE.20) THEN
667 + ((OAXIS(1)-PU(1+NPU))**2+PXMDPI(OAXIS(2)-PU(2+NPU))**2)
670 IF (COSVAL.GE.COSR)THEN
675 PNEW(MU) = PNEW(MU) + PP(MU+NPP)
679 + + PP(4+NPP)/(PP(4+NPP)+PNEW(4))*(PP(1+NPP)-PNEW(1))
681 + + PP(4+NPP)/(PP(4+NPP)+PNEW(4))
682 + *PXMDPI(PP(2+NPP)-PNEW(2))
683 PNEW(4)=PNEW(4)+PP(4+NPP)
689 *** If there are particles in the cone, calc new jet axis
694 NORMSQ = NORMSQ + PNEW(MU)**2
701 NAXIS(MU) = PNEW(MU)/NORM
707 *********************************************************************
708 *CMZ : 2.00/00 10/01/95 10.17.57 by P. Schleper
709 *CMZ : 1.06/00 28/02/94 15.44.44 by P. Schleper
713 SUBROUTINE PXUVEC(NTRAK,PP,PU,IERR)
715 *** Routine to calculate unit vectors PU of all particles PP
718 PARAMETER (MXTRAK=5000)
720 DOUBLE PRECISION PP(4,MXTRAK)
721 DOUBLE PRECISION PU(3,MXTRAK)
731 WRITE(6,*)' PXCONE: An input particle has zero mod(p)'
736 PU(MU,N)=PP(MU,N)/MAG
741 *CMZ : 1.06/00 28/02/94 15.44.44 by P. Schleper
743 C-----------------------------------------------------------------------
744 SUBROUTINE PXZERI(N,A)
750 *CMZ : 1.06/00 28/02/94 15.44.44 by P. Schleper
752 C-----------------------------------------------------------------------
753 C This is a set of routines written by Mike Seymour to provide the
754 C services presumably normally provided by standard OPAL routines
755 C PXZERV zeroes a vector
756 C PXZERI zeroes a vector of integers
757 C PXNORV normalizes a vector
758 C PXADDV adds two vectors
759 C PXSORV sorts a vector (copied from HERWIG)
760 C PXANG3 finds the angle (and its cosine) between two vectors
761 C PXMDPI moves its argument onto the range [-pi,pi)
762 C-----------------------------------------------------------------------
763 SUBROUTINE PXZERV(N,A)
765 DOUBLE PRECISION A(N)
771 C-----------------------------------------------------------------------
772 SUBROUTINE PXADDV(N,A,B,C,ITERR)
774 DOUBLE PRECISION A(N),B(N),C(N)
779 *CMZ : 1.06/00 28/02/94 15.44.44 by P. Schleper
781 C-----------------------------------------------------------------------
782 SUBROUTINE PXANG3(A,B,COST,THET,ITERR)
784 DOUBLE PRECISION A(3),B(3),C,COST,THET
785 C=(A(1)**2+A(2)**2+A(3)**2)*(B(1)**2+B(2)**2+B(3)**2)
788 COST=(A(1)*B(1)+A(2)*B(2)+A(3)*B(3))*C
791 *CMZ : 1.06/00 14/03/94 15.41.57 by P. Schleper
792 *-- Author : P. Schleper 28/02/94
793 LOGICAL FUNCTION PXNEW(TSTLIS,JETLIS,NTRAK,NJET)
795 INTEGER MXTRAK,MXPROT
796 PARAMETER (MXTRAK=5000,MXPROT=5000)
798 *** Note that although JETLIS is assumed to be a 2d array, it
799 *** it is used as 1d in this routine for efficiency
800 LOGICAL TSTLIS(MXTRAK),JETLIS(MXPROT*MXTRAK)
801 *** Checks to see if TSTLIS entries correspond to a jet already found
802 *** and entered in JETLIS
812 IF(TSTLIS(N).NEQV.JETLIS(IN)) THEN
824 *CMZ : 1.06/00 14/03/94 15.41.57 by P. Schleper
825 *-- Author : P. Schleper 28/02/94
826 LOGICAL FUNCTION PXSAME(LIST1,LIST2,N)
828 LOGICAL LIST1(*),LIST2(*)
830 *** Returns T if the first N elements of LIST1 are the same as the
831 *** first N elements of LIST2.
836 IF ( LIST1(I).NEQV.LIST2(I) ) THEN
843 *CMZ : 1.06/00 28/02/94 15.44.44 by P. Schleper
845 C-----------------------------------------------------------------------
848 C---RETURNS PHI, MOVED ONTO THE RANGE [-PI,PI)
849 DOUBLE PRECISION PXMDPI,PHI,PI,TWOPI,THRPI,EPS
850 PARAMETER (PI=3.141592654,TWOPI=6.283185307,THRPI=9.424777961)
851 PARAMETER (EPS=1E-15)
853 IF (PXMDPI.LE.PI) THEN
854 IF (PXMDPI.GT.-PI) THEN
856 ELSEIF (PXMDPI.GT.-THRPI) THEN
859 PXMDPI=-MOD(PI-PXMDPI,TWOPI)+PI
861 ELSEIF (PXMDPI.LE.THRPI) THEN
864 PXMDPI=MOD(PI+PXMDPI,TWOPI)-PI
866 100 IF (ABS(PXMDPI).LT.EPS) PXMDPI=0