*$ CREATE USRMED.FOR
*COPY USRMED
*                                                                      *
*=== usrmed ===========================================================*
*                                                                      *
      SUBROUTINE USRMED ( IJ, EKSCO, PLA, WEE, MREG, NEWREG, XX, YY, ZZ,
     &                    TXX, TYY, TZZ )

      INCLUDE '(DBLPRC)'
      INCLUDE '(DIMPAR)'
      INCLUDE '(IOUNIT)'
*
*----------------------------------------------------------------------*
*                                                                      *
*     Copyright (C) 1991-2005      by    Alfredo Ferrari & Paola Sala  *
*     All Rights Reserved.                                             *
*                                                                      *
*                                                                      *
*     USeR MEDium dependent directives:                                *
*                                                                      *
*     Created on  10  may  1996    by    Alfredo Ferrari & Paola Sala  *
*                                                   Infn - Milan       *
*                                                                      *
*     Last change on   29-may-96   by    Alfredo Ferrari               *
*                                                                      *
*     Input variables:                                                 *
*             ij = particle id                                         *
*          Eksco = particle kinetic energy (GeV)                       *
*            Pla = particle momentum (GeV/c)                           *
*            Wee = particle weight                                     *
*           Mreg = (original) region number                            *
*         Newreg = (final)    region number                            *
*       Xx,Yy,Zz = particle position                                   *
*    Txx,Tyy,Tzz = particle direction                                  *
*                                                                      *
*     The user is supposed to change only WEE if MREG = NEWREG and     *
*     WEE, NEWREG, TXX, TYY, TZZ if MREG .NE. NEWREG                   *
*                                                                      *
*----------------------------------------------------------------------*
* Written by Niels Bassler <bassler@phys.au.dk>
* This is an example of a RiFi modulator. Similar technique can 
* also be applied to a to ridge filters and range modulator wheels.
* In latter case, range sampling could be derived from a
* random number distribution.
* This is also described in http://dx.doi.org/10.1016/j.radonc.2010.02.026
* (Please cite, when used. :)
* Note, the data below is _NOT_ the GSI RiFi.
*
      REAL DZ, YYY
      INTEGER NNN

      REAL RIFI_Z0
      SAVE RIFI_Z0

      LOGICAL LFIRST
      SAVE LFIRST
      DATA LFIRST /.TRUE./

      LOGICAL LACT
      SAVE LACT
      DATA LACT /.TRUE./


      REAL DZZ(47)
      SAVE DZZ
      REAL YB(47)
      SAVE YB

* RIFI data (all in mm!) 
      DATA YB /0.0000000,0.0000015,0.0000024,0.0000043,0.0000053, 
     $     0.0003455,0.0003909,0.0004102,0.0004102,0.0005862,0.0016394,
     $     0.0078486,0.0273136,0.0555144,0.0879239,0.1096194,0.1370153,
     $     0.1568140,0.1858080,0.2057660,0.2366269,0.2570168,0.2932704,
     $     0.3116073,0.3501200,0.3696446,0.4078143,0.4290243,0.4660987,
     $     0.4861284,0.5223476,0.5412143,0.5734729,0.5906643,0.6186960,
     $     0.6332235,0.6548406,0.6667434,0.6788370,0.6917452,0.7069458,
     $     0.7199101,0.7391596,0.7472127,0.7472127,0.7481090,0.7500000 / 

      
      DATA DZZ /0.00,0.08,0.17,0.25,0.34,0.42,0.50,0.59,0.67,0.76,0.84,
     $     0.92,1.01,1.09,1.18,1.26,1.34,1.43,1.51,1.60,1.68,1.76,1.85,
     $     1.93,2.02,2.10,2.18,2.27,2.35,2.44,2.52,2.61,2.69,2.77,2.86,
     $     2.94,3.03,3.11,3.19,3.28,3.36,3.45,3.53,3.61,3.70,3.78,3.87 /

      


*     $FLUPRO/flutil/ldpm3qmd usrmed_rifi_new3mm.f -o flukadpm3_rifi_new3mm
* return message from first call
      IF (LFIRST) THEN
         WRITE(LUNOUT,*) '#NB USRMED_RIFI_PHSP 3mm RIFI ENABLED'
         WRITE(LUNOUT,*) 'NB IMPORTANT: CHECK ALSO FOR RIFI APPLIED!' 
         LFIRST = .FALSE.
      ENDIF
      
****************************************************
*     Set start position of RIFI here.
      RIFI_Z0 = -65.3
*****************************************************
      
      
* -----------------------------------------------------
*     Period of ripples is 0.15 cm 
*     
      
      IF (MREG .NE. NEWREG) THEN
         IF ((ZZ .GE. RIFI_Z0 - 0.00001) .AND. 
     &        (ZZ .LE. RIFI_Z0 + 0.00001) .AND.
     &        (TZZ .GT. 0)) THEN
      WRITE(LUNOUT,*) '_NB ',ZZ,RIFI_Z0-0.00001,RIFI_Z0+0.00001
            IF (LACT) THEN
               WRITE(LUNOUT,*) 'NB USRMED_RIFI 3mm RIFI APPLIED at:',ZZ
               LACT = .FALSE.
            END IF
            
            YYY = ABS(YY - (0.15 * AINT(YY/0.15)))
            IF (YYY .GT. 0.075) THEN
               YYY = 0.15 - YYY
            ENDIF
            
            NNN = 1
 10         IF (YYY .GT. YB(NNN)*0.1) THEN
*      WRITE(LUNOUT,*) '_NB_ LOOP ', NNN, YYY, YB(NNN)
               NNN = NNN + 1
               GOTO 10
            ENDIF

* secure array index, should never happen:
            IF (NNN .GT. 48) THEN
               WRITE(LUNOUT,*) '_NB_ Not good:', YY , YYY, NNN
               NNN = 48 
            ENDIF
               

* At YB = 0.75 mm do no shift, at YB = 0.00 shift max. 

            DZ = 0.387 - (DZZ(NNN)*0.1)
*            WRITE(LUNOUT,*) '_NB_ ', YY , YYY, NNN, YB(NNN)
*            WRITE(LUNOUT,*) '_NB NNN, DZZ, DZ', NNN, DZZ(NNN), DZ 

            ZZ = ZZ + DZ
         ENDIF             
      ENDIF       

*      WRITE(LUNOUT,*) "_NB Done/Return."
                  
      RETURN
*=== End of subroutine Usrmed =========================================*      
      END
      
      
