#include "cfx5ext.h"
dllexport(pt_collision)
      SUBROUTINE PT_COLLISION (NLOC,NRET,NARG,RET,ARG,CRESLT,
     &                        CZ,DZ,IZ,LZ,RZ)
CC
CD User routine: template for particle-particle collision user routine
CC
CC --------------------
CC        Input
CC --------------------
CC
CC  NLOC   - number of entities
CC  NRET   - length of return stack
CC  NARG   - length of argument stack
CC  ARG    - argument values
CC
CC --------------------
CC      Modified
CC --------------------
CC
CC --------------------
CC        Output
CC --------------------
CC
CC  RET    - return values
CC
CC --------------------
CC       Details
CC --------------------
CC
CC======================================================================
C
C ------------------------------
C     Preprocessor includes
C ------------------------------
C
#include "cfd_sysdep.h"
#include "cfd_constants.h"
C
C ------------------------------
C        Argument list
C ------------------------------
C
      INTEGER NARG, NRET, NLOC
C
      CHARACTER*(4) CRESLT
C
      REAL ARG(NLOC,NARG), RET(NLOC,NRET)
C
      INTEGER IZ(*)
      CHARACTER CZ(*)*(1)
      DOUBLE PRECISION DZ(*)
      LOGICAL LZ(*)
      REAL RZ(*)
C
C=======================================================================
C
C ---------------------------
C    Executable Statements
C ---------------------------
C
C=======================================================================
C
C     Return variables
C     ---------------------------------
C
C     Particle diameter                                    : RET(1,1)
C     Particle velocity                                    : RET(1,2)
C
C     Argument variables
C     -----------------------------------
C
C     Particle density                                     : ARG(1,1)
C     Particle position                                    : ARG(1,2)
C
C     Particle diameter                                    : ARG(1,5)
C     Particle mean diameter                               : ARG(1,6)
C     Particle RMS diameter                                : ARG(1,7)
C
C     Particle velocity                                    : ARG(1,8)
C     Particle averaged velocity                           : ARG(1,11)
C     Particle RMS velocity                                : ARG(1,14)
C
C     Particle number density                              : ARG(1,17)
C     Particle integration time step                       : ARG(1,18)
C     Particle turbulent Stokes number                     : ARG(1,19)
C
C=======================================================================
C
C-----------------------------------------------------------------------
C     Calculate particle-particle collision according to Sommerfeld
C     model
C-----------------------------------------------------------------------
C
      CALL COLL_SOMMERFELD (RET(1,1),RET(1,2),
     &                      ARG(1,1),ARG(1,2),
     &                      ARG(1,5),ARG(1,6),ARG(1,7),
     &                      ARG(1,8),ARG(1,11),ARG(1,14),
     &                      ARG(1,17),ARG(1,18),ARG(1,19))
C
      END
C
C
C===============================================================
C===============================================================
C    User routine follows
C===============================================================
C===============================================================
C
C
      SUBROUTINE COLL_SOMMERFELD(dP1_prime, v1_prime,
     &                           rhoP, P1globcoords,
     &                           dP1, dPmean, sigma_dP,
     &                           v1, v1_mean, sigma_P,
     &                           nP, Delta_t, St_t)
CC
CC Routine for the stochastic particle-particle collision model
CC by B. Oesterle and M. Sommerfeld
CC
CC Implementation by ANSYS CFX Germany
CC and University of the Federal Armed Forces Munich
CC
C
      IMPLICIT NONE
C 
C ---------------------------
C    Preprocessor includes
C ---------------------------
C
#include "cfd_constants.h"
C
C ------------------------------
C        Argument list
C ------------------------------
C
C  dP1_prime:    diameter of particle P1 after collision
C  v1_prime:     post-collision velocity of P1 in global xyz-system
C
C ------------------------
C    Output arguments
C ------------------------
C
      REAL     dP1_prime, v1_prime(3)
C
C ------------------------
C     Input arguments
C ------------------------
C
CI rhoP:          density of particles
CI P1globcoords:  global coordinates of particle P1
CI dP1:           diameter of particle 1
CI dPmean:        local mean diameter of particles in previous 
CI                Lagrange-iteration
CI sigma_dP:      local standard deviation of particle diameters 
CI                in previous iteration
CI v1:            velocity of particle P1
CI v1_mean:       local mean velocity of particle P1 during the 
CI                previous  Lagr.-iteration
CI sigma_P:       local average fluctuation velocity
CI nP:            particle number density
CI Delta_t:       Lagrangian time step (has to be passed back 
CI                to the solver)
CI St_t:          turbulent Stokes number
 
      REAL    rhoP, P1globcoords(3), dP1, dPmean, sigma_dP, v1(3),
     &        v1_mean(3), sigma_P(3), nP, Delta_t, St_t
C
C -----------------------------------------------------
C     Model constants
C -----------------------------------------------------
C
CI e:      coefficient of restitution
CI mu_fs:  coefficient of static friction
CI mu_fk:  coefficient of kinetic friction
C
      REAL    e, mu_fs, mu_fk
C
C ------------------------------
C        External routines
C ------------------------------
C
      REAL     ran_coll
      EXTERNAL ran_coll
C
C ------------------------------
C        Local Parameters
C ------------------------------
C
C ------------------------------
C        Local variables
C ------------------------------
C
      INTEGER IDIM, NDIM, ISEED, ACTION, WHAT
C
      REAL    R_St_t, v1_fluc(3), v2_fluc(3), xi(1), v2(3),
     &        zeta(1), dP2, chi(1), fc, Pc, psi, n_prime(3), t_prime(3),
     &        k_prime(3), unit1(3), P2localcoords(3), P2globcoords(3),
     &        rotmat(3,3), n(3), t(3), k(3), ktemp(3), O1O2(3),
     &        normO1O2, vPrel(3), v1ntk(3), v2ntk(3),
     &        vR(3), m1, m2, v1ntk_prime(3),
     &        vPrelAbs, ktempNorm, SN, FACT
C
C ---------------------------
C     Data Statements
C ---------------------------
C
C-----------------------------------------------------------------------
C     Initialisation of variables
C     --> Particle rotation is set to zero
C-----------------------------------------------------------------------
C
C---- Get ISEED from particle tracker
C       ACTION = 2 -> __pt_uf_save_data__
C       WHAT   = 1 -> __pt_uf_iseed__
C
      ACTION = 2
      WHAT   = 1
      CALL USER_PARTICLE_INFO(ACTION,WHAT,ISEED)
C
C---- Initialization
C
      e    = 0.95
      NDIM = 3
      SN   = 1.0E-20
C
C---- This yields a frictionless sliding collision:
C
      mu_fs = 0.0
      mu_fk = 0.0
C
C---- If number density is too small, there is no collision partner,
C     so return
C
      IF (nP .LT. 1.0e-6) THEN
         DO IDIM=1,NDIM
            v1_prime(IDIM)     = v1(IDIM)
         ENDDO
         dP1_prime = dP1
         Delta_t   = 1.0E+30
         GOTO 999
      ENDIF         
C
C-----------------------------------------------------------------------
C     Calculating the instantaneous velocity of the fictitious 
C     collision partner P2
C-----------------------------------------------------------------------
C
C---- Velocity correlation function
C
      R_St_t = EXP(-0.55 * St_t**0.4)
C
C---- Fluctuation velocity of P1 (assume NDIM = 3)
C
      CALL SET_A_BminusC(v1_fluc,v1,v1_mean,NDIM)
C
C---- Fluctuation velocity of the virtual collision partner P2
C
      DO IDIM = 1,NDIM
        CALL GET_RANDOM_NORM(XI,IONE,ZERO,ONE,-FIVE,FIVE,ISEED)
        v2_fluc(IDIM) = R_St_t*v1_fluc(IDIM)
     &                + sigma_P(IDIM)*SQRT(1-R_St_t**2 )*xi(1)
      ENDDO
C
C---- Instantaneous velocity of fictitious particle P2
C
      CALL SET_A_BplusC(v2,v1_mean,v2_fluc,NDIM)
C
C-----------------------------------------------------------------------
C     Collision frequency and collision probability
C-----------------------------------------------------------------------
C
C---- Calculation of diameter of fictitious particle P2
C
      CALL GET_RANDOM_NORM(CHI,IONE,ZERO,ONE,-FIVE,FIVE,ISEED)
      dP2 = dPmean + sigma_dP * chi(1)
C
C---- Collision frequency
C
      CALL SET_A_BminusC(vPrel,v1,v2,NDIM)
      CALL SET_R_normA(vPrelAbs,vPrel,NDIM)
C
      fc = PI/4. * (dP1 + dP2)**2 * vPrelAbs * nP
C
C---- Collision probability (using the "old" time step)
C
      Pc = 1 - EXP(-fc * Delta_t)
C
C----  Check if Lagrangian time step exceeds the stability limit.
C      --> Adjust Lagrangian time step, but do not let it go too small
C
      IF (FC .GT. SN .AND. Delta_t .GT. 0.05/(fc + SN)) THEN
        Delta_t = 0.04/fc
        Delta_t = MAX(Delta_t,1.0E-10)
      ELSE
        Delta_t = 1.0E+30
      ENDIF
C
C---- Decide whether a collision takes place or not
C
      psi = ran_coll(ISEED)
C
C---- Psi > Collision probability:
C     --> No collision takes place
C
      IF (psi .GT. Pc) THEN
        DO IDIM=1,NDIM
          v1_prime(IDIM) = v1(IDIM)
        ENDDO
        dP1_prime = dP1
        GOTO 999
      ENDIF
C
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C     Collision takes place
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C
C---- Compute location of the virtual collision partner
C     --> Calculation of the local coordinate system base vectors
C
      CALL SET_A_BdivbyR(n_prime,vPrel,vPrelAbs,NDIM)

      IF (n_prime(1).EQ.1.0 .AND. 
     &    n_prime(2).EQ.0.0 .AND.
     &    n_prime(3).EQ.0.0) THEN
        t_prime(1) = ZERO
        t_prime(2) = ONE
        t_prime(3) = ZERO
C
        k_prime(1) = ZERO
        k_prime(2) = ZERO
        k_prime(3) = ONE
      ELSE
        unit1(1)   = ONE
        unit1(2)   = ZERO
        unit1(3)   = ZERO
C
        CALL VCROSS(t_prime, n_prime, unit1)
        CALL VCROSS(k_prime, n_prime, t_prime)
      ENDIF
C
C---- Calculation of the position of particle P2 in local coordinates
C
      P2localcoords(2) = (dP1 + dP2)*(ran_coll(iseed) - HALF)
      P2localcoords(3) = (dP1 + dP2)*(ran_coll(iseed) - HALF)
C
      DO WHILE (P2localcoords(2)**2 + P2localcoords(3)**2 .GT.
     &          QUARTER*(dP1 + dP2)**2)
        P2localcoords(2) = (dP1 + dP2)*(ran_coll(iseed) - HALF)
        P2localcoords(3) = (dP1 + dP2)*(ran_coll(iseed) - HALF)
      ENDDO
C
      FACT = QUARTER*(dP1 + dP2)**2 - P2localcoords(2)**2 
     &                              - P2localcoords(3)**2
C
      P2localcoords(1) = SQRT(MAX(0.0,FACT))
C
C---- Inverse transformation back to the global coordinate system
C     comprising a rotation and a parallel translation by the
C     position vector of particle P1

C---- Rotation matrix
C
      DO IDIM=1,NDIM
        rotmat(IDIM,1) = n_prime(IDIM)
        rotmat(IDIM,2) = t_prime(IDIM)
        rotmat(IDIM,3) = k_prime(IDIM)
      ENDDO
C
C---- Rotation and parallel translation
C
      CALL MatrProd(rotmat,P2localcoords,P2globcoords) 
C
      CALL SET_A_AplusB(P2globcoords,P1globcoords,NDIM)
C
C-----------------------------------------------------------------------
C     Deterministic calculation of a binary collision
C-----------------------------------------------------------------------
C
C---- Creating the new local coordinate system:
C
      CALL SET_A_BminusC(O1O2,P2globcoords,P1globcoords,NDIM)
C
      CALL SET_R_normA(normO1O2,O1O2,NDIM)
      CALL SET_A_BdivbyR(n,O1O2,normO1O2,NDIM)
C
      CALL VCROSS(ktemp, n, vPrel)
C
      CALL SET_R_normA(ktempNorm,ktemp,NDIM)
      CALL SET_A_BdivbyR(k,ktemp,ktempNorm,NDIM)
C
      CALL VCROSS(t, k, n) 
C
C---- Forward transformation of all velocities to the local coordinate
C     system comprising (only) a rotation between global and local 
C     coordinate axes
C
      DO IDIM=1,NDIM
        rotmat(1,IDIM) = n(IDIM)
        rotmat(2,IDIM) = t(IDIM)
        rotmat(3,IDIM) = k(IDIM)
      ENDDO
C
C---- Rotation of all 4 velocities
C
      CALL MatrProd(rotmat,v1,v1ntk) 
      CALL MatrProd(rotmat,v2,v2ntk) 
C
C---- Calculation of relative velocity vR in contact point M
C
      vR(1) = v1ntk(1) - v2ntk(1)
      vR(2) = v1ntk(2) - v2ntk(2) 
      vR(3) = v1ntk(3) - v2ntk(3) 
C
C---- Calculation of mass of particle P1 and P2, respectively
C
      m1 = rhoP * 4./3. * PI * dP1**3
      m2 = rhoP * 4./3. * PI * dP2**3
C
C---- Case destinction between non-sliding and sliding collision
C
      IF (SQRT(vR(2)**2+vR(3)**2).LE. 7./2.*(1+e)*mu_fs*ABS(vR(1))) THEN
C
C---- NON-sliding collision
C
         v1ntk_prime(1) = v1ntk(1) - m2/(m1 + m2) * (1 + e) * vR(1)
         v1ntk_prime(2) = v1ntk(2) - 2./7. * m2/(m1 + m2) * vR(2)
         v1ntk_prime(3) = v1ntk(3) - 2./7. * m2/(m1 + m2) * vR(3)
      ELSE
C
C---- SLIDING collision occurs
C
         v1ntk_prime(1) = v1ntk(1) - m2/(m1 + m2)*(1 + e)*vR(1)
         v1ntk_prime(2) = v1ntk(2) - mu_fk * vR(2)/(SQRT(vR(2)**2 
     &                  + vR(3)**2))*m2/(m1 + m2)*(1 + e)*ABS(vR(1))
         v1ntk_prime(3) = v1ntk(3) - mu_fk * vR(3)/(SQRT(vR(2)**2 
     &                  + vR(3)**2))*m2/(m1 + m2)*(1 + e)*ABS(vR(1))
      ENDIF
C
C---- Inverse transformation back to the global coordinate system
C     comprising only a rotation
C
      DO IDIM=1,NDIM
        rotmat(IDIM,1) = n(IDIM)
        rotmat(IDIM,2) = t(IDIM)
        rotmat(IDIM,3) = k(IDIM)
      ENDDO
C
C---- Rotation of velocity
C
      CALL MatrProd(rotmat,v1ntk_prime,v1_prime)
C
C---- Particle diameter is not changed by this routine
C
      dP1_prime = dP1
C
C -----------------------------------------------------------
C ----  dP1_prime and v1_prime are the results of the    ----
C ----  routine and are passed back to the solver.       ----
C -----------------------------------------------------------
C
 999  CONTINUE
C
C---- Store modified ISEED and return to particle solver
C       ACTION = 1 -> __pt_uf_get_data__
C       WHAT   = 1 -> __pt_uf_iseed__
C
      ACTION = 1
      WHAT   = 1
      CALL USER_PARTICLE_INFO(ACTION,WHAT,ISEED)
C
      RETURN
      END
C
C=======================================================================
C=======================================================================
C=======================================================================
C
      FUNCTION RAN_COLL(ISEED)
CC
CD    Interface to random number generator (GET_RANDOM)
CC
C ------------------------------
C        Argument list
C ------------------------------
C
      INTEGER ISEED
      REAL    RAN_COLL
C
C ------------------------------
C        Local variables
C ------------------------------
C
      REAL    RAND(1)
C
      CALL GET_RANDOM(RAND,1,ISEED)
C
      RAN_COLL = RAND(1)
      END
