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

      SUBROUTINE UPHI(IENTRY,LVL)
C                                VERSION 4.00  --  26 JAN 1986/1900
C******************************************************************
C   UPHI STANDS FOR 'UNIFORM PHI DISTRIBUTION'.
C   SET COORDINATES FOR NEW PARTICLE OR RESET DIRECTION COSINES OF
C   OLD ONE.  GENERATE RANDOM AZIMUTH SELECTION AND REPLACE THE
C   DIRECTION COSINES WITH THEIR NEW VALUES.
C******************************************************************
*KEEP,EPCONT.
      COMMON/EPCONT/   EDEP,RATIO,TSTEP,TUSTEP,USTEP,TVSTEP,VSTEP,IDISC,
     *                 IROLD,IRNEW,RHOFAC, EOLD,ENEW,EKE,ELKE,BETA2,GLE,
     *                 TSCAT,IAUSFL
      DOUBLE PRECISION EDEP,RATIO
      REAL             TSTEP,TUSTEP,USTEP,TVSTEP,VSTEP,RHOFAC,EOLD,ENEW,
     *                 EKE,ELKE,BETA2,GLE,TSCAT
      INTEGER          IDISC,IROLD,IRNEW,IAUSFL(29)
*KEEP,RANDPA.
      COMMON /RANDPA/  FAC,U1,U2,RD,NSEQ,ISEED,KNOR
      DOUBLE PRECISION FAC,U1,U2
      REAL             RD(3000)
      INTEGER          ISEED(103,10),NSEQ
      LOGICAL          KNOR
*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
*KEEP,STACKE.
      COMMON/STACKE/   E,TIME,X,Y,Z,U,V,W,DNEAR,IQ,IGEN,IR,IOBS,LPCTE,NP
      DOUBLE PRECISION E(60),TIME(60)
      REAL             X(60),Y(60),Z(60),U(60),V(60),W(60),DNEAR(60)
      INTEGER          IQ(60),IGEN(60),IR(60),IOBS(60),LPCTE(60),NP
*KEND.
      COMMON/UPHIIN/SINC0,SINC1,SIN0(20002),SIN1(20002)
      COMMON/UPHIOT/THETA,SINTHE,COSTHE,SINPHI, COSPHI,PI,TWOPI,PI5D2
      SAVE A,B,C
      IF((IENTRY.EQ.2))GO TO1070
      IF((IENTRY.EQ.3))GO TO1080
1090  LTHETA=SINC1*THETA+SINC0
      SINTHE=SIN1(LTHETA)*THETA+SIN0(LTHETA)
      CTHET=PI5D2-THETA
      LCTHET=SINC1*CTHET+SINC0
      COSTHE=SIN1(LCTHET)*CTHET+SIN0(LCTHET)
C   USE THE FOLLOWING ENTRY IF SINTHE AND COSTHE ARE ALREADY KNOWN.
C   SELECT PHI UNIFORMLY OVER THE INTERVAL (0,TWO PI). THEN USE
C   PWLF OF SIN FUNCTION TO GET SIN(PHI) AND COS(PHI).  THE COSINE
C   IS GOTTEN BY COS(PHI)=SIN(9*PI/4 - PHI).
1070  CALL RMMAR(RNNO38,1,2)
      PHI=RNNO38*TWOPI
      LPHI=SINC1*PHI+SINC0
      SINPHI=SIN1(LPHI)*PHI+SIN0(LPHI)
      CPHI=PI5D2-PHI
      LCPHI=SINC1*CPHI+SINC0
      COSPHI=SIN1(LCPHI)*CPHI+SIN0(LCPHI)
C   USE THE FOLLOWING ENTRY FOR THE SECOND OF TWO PARTICLES WHEN WE
C   KNOW TWO PARTICLES HAVE A RELATIONSHIP IN THEIR CORRECTIONS.
C   NOTE: SINTHE AND COSTHE CAN BE CHANGED OUTSIDE THROUGH COMMON.
C   LVL IS A PARAMETER TELLING WHICH PARTICLES TO WORK WITH.
C   THETA (SINTHE AND COSTHE) ARE ALWAYS RELATIVE TO THE DIRECTION
C   OF THE INCIDENT PARTICLE BEFORE ITS DIRECTION WAS ADJUSTED.
C   THUS WHEN TWO PARTICLES NEED TO HAVE THEIR DIRECTIONS COMPUTED,
C   THE ORIGINAL INCIDENT DIRECTION IS SAVED IN THE VARIABLE A,B,C
C   SO THAT IT CAN BE USED ON BOTH CALLS.
C   LVL=1 -- OLD PARTICLE, SAVE ITS DIRECTION AND ADJUST IT
C   LVL=2 -- NEW PARTICLE. ADJUST DIRECTION USING SAVED A,B,C
C   LVL=3 -- BREMSSTRAHLUNG GAMMA.  SAVE ELECTRON DIRECTION (NEXT
C   TO TOP OF STACK), AND THEN ADJUST GAMMA DIRECTION.
1080  IF (LVL.EQ.2) GO TO1100
      IF((LVL.EQ.3))GO TO1110
1120  A=U(NP)
      B=V(NP)
      C=W(NP)
      GO TO 1130
1110  A=U(NP-1)
      B=V(NP-1)
      C=W(NP-1)
1100  X(NP)=X(NP-1)
      Y(NP)=Y(NP-1)
      Z(NP)=Z(NP-1)
      LPCTE(NP)=LPCTE(NP-1)
      IR(NP)=IR(NP-1)
      DNEAR(NP)=DNEAR(NP-1)
      TIME(NP)=TIME(NP-1)
      IGEN(NP)=IGEN(NP-1)
      IOBS(NP)=IOBS(NP-1)
1130  SINPS2=A*A+B*B
      IF (SINPS2.LT.1.0E-10) THEN
       U(NP)=SINTHE*COSPHI
       V(NP)=SINTHE*SINPHI
       W(NP)=C*COSTHE
      ELSE
       SINPSI=SQRT(SINPS2)
       US=SINTHE*COSPHI
       VS=SINTHE*SINPHI
       SINDEL=B*(1./SINPSI)
       COSDEL=A*(1./SINPSI)
       U(NP)=C*COSDEL*US-SINDEL*VS+A*COSTHE
       V(NP)=C*SINDEL*US+COSDEL*VS+B*COSTHE
       W(NP)=-SINPSI*US+C*COSTHE
      END IF
      RADINV=1.5-0.5*(U(NP)**2+V(NP)**2+W(NP)**2)
      U(NP)=U(NP)*RADINV
      V(NP)=V(NP)*RADINV
      W(NP)=W(NP)*RADINV
      RETURN
      END