source: trunk/MagicSoft/Simulation/Corsika/Mmcs/jadach.f

      SUBROUTINE JADACH( ECMJAD,JADFLG )
 
C-----------------------------------------------------------------------
C  JADACH (FILTER)
C
C  ADJUSTS THE RAPIDITIES OF ALL SECONDARIES SUCH THAT
C  ENERGY AND LONGITUDINAL MOMENTUM ARE CONSERVED AT THE SAME TIME
C  THE ALGORITHM IS TAKEN FROM S.JADACH, COM.PHYS.COMM. 9 (1975) 297
C  THE ROUTINE MUST BE CALLED AFTER THE PT IS CONSERVED AND BEFORE
C  THE TRANSFORMATION TO THE LAB SYSTEM IS DONE
C  THIS SUBROUTINE IS CALLED FROM HDPM
C  ARGUMENTS:
C   ECMJAD = CM ENERGY IN THE PROJECTILE -- GNU*NUCLEONS SYSTEM
C   JADFLG = 0  JADACH FILTER CORRECTLY ENDED
C          = 1  BAD RAPIDITIES, SELECT RAPIDITIES AGAIN
C          =-1  SUM OF TRANSVERSE MASSES EXCEEDS AVAILABLE CM ENERGY
C-----------------------------------------------------------------------
 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
*KEEP,INTER.
      COMMON /INTER/   AVCH,AVCH3,DC0,DLOG,DMLOG,ECMDIF,ECMDPM,ELAB,
     *                 FNEUT,FNEUT2,GNU,PLAB,POSC2,POSC3,POSN2,POSN3,
     *                 RC3TO2,S,SEUGF,SEUGP,SLOG,SLOGSQ,SMLOG,
     *                 WIDC2,WIDC3,WIDN2,WIDN3,YCM,YY0,ZN,
     *                 IDIF,ITAR
      DOUBLE PRECISION AVCH,AVCH3,DC0,DLOG,DMLOG,ECMDIF,ECMDPM,ELAB,
     *                 FNEUT,FNEUT2,GNU,PLAB,POSC2,POSC3,POSN2,POSN3,
     *                 RC3TO2,S,SEUGF,SEUGP,SLOG,SLOGSQ,SMLOG,
     *                 WIDC2,WIDC3,WIDN2,WIDN3,YCM,YY0,ZN
      INTEGER          IDIF,ITAR
*KEEP,NEWPAR.
      COMMON /NEWPAR/  EA,PT2,PX,PY,TMAS,YR,ITYP,
     *                 IA1,IA2,IB1,IB2,IC1,IC2,ID1,ID2,IE1,IE2,IF1,IF2,
     *                 IG1,IG2,IH1,IH2,II1,II2,IJ1,NTOT
      DOUBLE PRECISION EA(3000),PT2(3000),PX(3000),PY(3000),TMAS(3000),
     *                 YR(3000)
      INTEGER          ITYP(3000),
     *                 IA1,IA2,IB1,IB2,IC1,IC2,ID1,ID2,IE1,IE2,IF1,IF2,
     *                 IG1,IG2,IH1,IH2,II1,II2,IJ1,NTOT
*KEEP,PAM.
      COMMON /PAM/     PAMA,SIGNUM
      DOUBLE PRECISION PAMA(6000),SIGNUM(6000)
*KEEP,RUNPAR.
      COMMON /RUNPAR/  FIXHEI,THICK0,HILOECM,HILOELB,
     *                 STEPFC,NRRUN,NSHOW,PATAPE,MONIIN,
     *                 MONIOU,MDEBUG,NUCNUC,
     *                 CETAPE,
     *                 SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
     *                 N1STTR,MDBASE,
     *                 DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
     *                 FIX1I,FMUADD,FNKG,FPRINT,FDBASE
     *                ,GHEISH,GHESIG
      COMMON /RUNPAC/  DSN,HOST,USER
      DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB
      REAL             STEPFC
      INTEGER          NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC,
     *                 SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
     *                 N1STTR,MDBASE
      INTEGER          CETAPE
      CHARACTER*79     DSN
      CHARACTER*20     HOST,USER
 
      LOGICAL          DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
     *                 FIX1I,FMUADD,FNKG,FPRINT,FDBASE
     *                ,GHEISH,GHESIG
*KEND.
 
      DIMENSION YRJAD(3000)
      DATA EPS / 1.D-7 /
C-----------------------------------------------------------------------
 
      IF ( DEBUG ) WRITE(MDEBUG,*) 'JADACH: NTOT=',NTOT
 
      JADFLG = 0
 
C  SUM UP TRANSVERSE MOMENTA AND COMPARE WITH AVAILABLE C.M.ENERGY
      STMAS = 0.D0
      ECMI  = 1.D0 / ECMJAD
      DO  4  I = 1,NTOT
        STMAS    = STMAS + TMAS(I)
        YRJAD(I) = YR(I)
  4   CONTINUE
      REST  = ( ECMJAD - STMAS ) * ECMI
      IF ( REST .LE. 0.D0 ) THEN
C  SUMMED TRANSVERSE MASS > AVAILABLE C.M. ENERGY
        JADFLG = -1
        RETURN
      ENDIF
      FACT   = 1.5D0 / REST
      AA     = 1.D0
      DIFOLD = 0.D0
      ICOUNT = 0
C  OPTIMIZATION LOOP TO DEFINE PARAMETER AA
  1   CONTINUE
      ICOUNT = ICOUNT + 1
      IF ( ICOUNT .GE. 50 ) GOTO 999
C  FORM SUMS S1 AND S2
      S1     = 0.D0
      S2     = 0.D0
      DO  5  I = 1,NTOT
        EXPO = EXP( AA * YR(I) )
        S1   = S1 + TMAS(I) * ECMI * EXPO
        S2   = S2 + TMAS(I) * ECMI / EXPO
  5   CONTINUE
      DIFF   = 0.1D0 * LOG(S1*S2)
C  ACCELERATING OF CONVERGENCE IF NO CHANGE OF SIGN IN DIFF
      IF ( DIFOLD*DIFF .GE. 0.D0 ) DIFF =  DIFF * FACT
      DIFOLD = DIFF
      IF ( DEBUG ) WRITE(MDEBUG,*) '   DIFF=',SNGL(DIFF)
      AA     = AA * MAX( 0.1D0, (1.D0 - DIFF) )
      IF ( ABS(DIFF) .GT. EPS ) GOTO 1
 
C  ITERATION HAS CONVERGED, CALCULATE PARAMETER BB
      BB     = 0.5D0 * LOG(S2/S1)
 
      IF ( DEBUG ) WRITE (MDEBUG,110) ICOUNT,STMAS,REST
  110 FORMAT('   ICOUNT, STMAS, REST = ',I5,2E13.5,/
     *       '   NUM   ITYP    TMAS          YR(OLD)         YR(NEW)')
C  CORRECT RAPIDITIES
      DO 10  I = 1,NTOT
        YR(I) = AA * YR(I) + BB
        IF ( DEBUG ) WRITE(MDEBUG,111) I,ITYP(I),TMAS(I),YRJAD(I),YR(I)
 111    FORMAT('  ',I4,I6,F12.5,2F16.8)
C  IMPOSSIBLE RAPIDITY, DETERMINE AGAIN THE RAPIDITIES
        IF ( ABS(YR(I)) .GT. LOG(ECMJAD/PAMA(ITYP(I))) ) GOTO 999
 10   CONTINUE
      RETURN
 
 
C  ERROR EXIT
 999  JADFLG = 1
C  NO CONVERGENCE AFTER 50 ITERATIONS OR IMPOSSIBLE RAPIDITY
      RETURN
 
      END