C This program calculates the mobility and temperature of ions
C in a weakly ionized gas as functions of E/n and T the neutral gas
C temperature by means of the approach described in: P. G. C. Almeida,
C M. S. Benilov, and G. V. Naidis, Calculation of Ion Mobilities by 
C Means of the Two-Temperature Displaced-Distribution Theory, J. Phys. 
C D: Appl. Phys. 35, No. 13, pp. 1577-1584 (2002). (This paper can be
C downloaded from the site of J. Phys. D: Appl. Phys. or from
C http://fisica.uma.pt/ingles/pessoal/Mikhail_Benilov/PDF/mobility.html)

C This source code is made available for downloading as a supplement
C to the executable file stored in the same directory. If you need
C to calculate a system ion-neutral gas that is not implemented
C in the executable, you can take this source code and add to it
C the corresponding energy-dependent cross section for momentum
C transfer in collisions of ions with neutral particles, QIA1, and
C masses of the corresponding ions and atoms. Please keep in mind
C that we have not tried to maintain the database in this source code
C up-to-date, so the database in the executable file includes a number
C of systems absent in this source code.

C The present version of the program was written by Pedro Almeida,
C pedroa@uma.pt

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON AA,BB,ama,ami
      parameter (ak=1.380658d-23,qe=1.60217733d-19,pi=3.141592654)
      INTEGER points
      CHARACTER*2 GAS
      CHARACTER*3 ions
      CHARACTER*1 fname
      CHARACTER*50 filename
      CHARACTER*3 escale

   31 WRITE (6,*)'Choose the gas (the following gases are available: He,
     N Ne, Ar, Kr, Xe, O2, N2): '
      READ *, GAS
      IF (GAS.NE.'Ar'.AND.GAS.NE.'He'.AND.GAS.NE.'Xe'
     ,   .AND.GAS.NE.'Ne'.AND.GAS.NE.'Kr'.and.gas.ne.'O2'.and.gas.ne.'N2
     '') T
     hHEN

         WRITE (6,*) 'The present version of the program is able to calc
     uulate the ion drift velocity, mobility, and temperature only for t
     hhe following gases: He, Ne, Ar, Kr, Xe, O2, N2 while you have spec
     iified ', GAS
         PRINT *, ' '
         GO TO 31
      END If

      IF (gas.eq.'He') THEN
   52     write (6,*) 'The program will make calculations for He+.'
          ions='He+'
      END IF

      IF (gas.eq.'Ne') THEN
   64     write (6,*) 'The program will make calculations for Ne+.'
          ions='Ne+'
      END IF

      IF (gas.eq.'Ar') THEN
   76     write (6,*) 'The program will make calculations for Ar+.'
          ions='Ar+'
      END IF

      IF (gas.eq.'Kr') THEN
   88     write (6,*) 'The program will make calculations for Kr+.'
          ions='Kr+'
      END IF

      IF (gas.eq.'Xe') THEN
          write (6,*) 'The program will make calculations for Xe+.'
          ions='Xe+'
      END IF

      IF (gas.eq.'O2') THEN
  112     write (6,*) 'The program will make calculations for O2+.'
          ions='O2+'
      END IF

      IF (gas.eq.'N2') THEN
          write (6,*) 'The program will make calculations for N2+.'
          ions='N2+'
      END IF

      CALL GASIONID(GAS,ions)

      PRINT *, ' '
   34 WRITE (6,*)'Input the lower limit of the E/n range you wish to cal
     cculate in Td: '
      READ *, Elow
      WRITE (6,*)'Input the upper limit of the E/n range you wish to cal
     cculate in Td: '
      READ *, Ehigh
      IF (Elow.lt.TINY(Elow).OR.Ehigh.lt.TINY(Ehigh)) THEN
           WRITE (6,*) 'The values of E/n, must be non-zero.'
           GO TO 34
      END IF
      IF (Elow.gt.Ehigh) THEN
           write (6,*) 'The lower limit of E/n must be smaller than the
     uupper limit.'
           GO TO 34
      END IF
      IF (Ehigh-Elow.lt.TINY(Ehigh-Elow)) THEN
         points=1
         step=0
         Else
           WRITE (6,*)'Input the number of E/n points you wish to calcul
     aate: '
           READ *, points
           IF (points.eq.0.or.points.eq.1) THEN
              write (6,*)'The number of points should be greater that 1.
     I If you would like to calculate just one point, then chose equal l
     oower and upper limits for E/n.'
              GO TO 34
           END IF
      END IF

      IF (points.eq.1) THEN
         escale='lin'
         GO TO 171
      END IF

  131 WRITE (6,*)'Choose beteween linear (lin) and logarithmic (log) sca
     lle for E/n: '
      read *,escale
      if (escale.ne.'lin'.and.escale.ne.'log') then
         GO TO 131
      end if

  171 WRITE (6,*)'Input T in K: '
      READ *, T
      IF (T.lt.TINY(T)) THEN
         write (6,*)'The value of T must be non-zero.'
         GO TO 171
      END IF

  138 WRITE (6,*)'Do you wish to specify the output file name? (y/n) '
      read *, fname
      IF (fname.ne.'y'.and.fname.ne.'n') THEN
         GO TO 138
      END IF
      IF (fname.eq.'y') THEN
         write (6,*)'Specify the file name where the results will be sto
     rred. It must be inserted in the form filename.extension: '
         read *, filename
         Else
         filename='mobility.out'
      END IF


      OPEN (9,file = filename)
      WRITE (9,'(10A18)') 'E/n(Td)', 'mui*n(1/(m*V*s))', 'Ti(K)'

      write (6,*)'Program running... please wait.'

      Ehigh=Ehigh*1.D-21
      Elow=Elow*1.D-21

      IF (escale.eq.'log') THEN
         Ehigh=dlog(Ehigh)
         Elow=dlog(Elow)
         step=(Ehigh-Elow)/(points-1)
         Else
         step=(Ehigh-Elow)/(points-1)
      END IF

C The integration limits for the weighted cross-sections are:
  163 A=1.D-3
      B=15.D0

C AM is the Mach number
C AMOB is the mobility*n
C VI is the ion drift velocity
C TI is the ion fluid temperature
C T is the neutral gas temperature

      z=ami*ama/(ami+ama) !z is the reduced mass
      AMOB0=3.*qe*pi/(16.*ak*T)*SQRT(2.*ak*T/(pi*z))/QIA11(A,B,T) !this is the first aproximation to the ion mobility
      m=0
      E=Elow
      DO WHILE (m.le.points-1)
         n=0 !n is the iteration counter
         IF (escale.eq.'log') THEN
            E=DEXP(E)
            else
            E=E
         END IF
         VI0=AMOB0*E
         TI0=T+ama/(3.*ak)*VI0**2
         X=100.D0
         DO WHILE (X.GT.1.D-6)
            Tia=(ama*Ti0+ami*T)/(ami+ama)
	    AM=SQRT(Z/(2.*ak*Tia))*VI0
            eta=SQRT(Tia/Ti0)*QIAM(AM,A,B,Tia)
            chi=Qiae(am,a,b,Tia)/Qiam(am,a,b,Tia)-2./3.*AM**2.
            AMOB=3.*qe*pi/(16.*ak*Ti0)*SQRT(2.*ak*Ti0/(pi*z))/eta
            VI1=AMOB*E
            TI1=T+ama/(3.*ak*CHI)*VI0**2
            X=ABS((VI1-VI0)/VI1)
            VI0=0.2*VI0+(1.D0-0.2)*VI1
            TI0=0.2*TI0+(1.D0-0.2)*TI1
            n=n+1
            Amob0=Amob
         END DO
         WRITE (9,'(10ES18.4)') E/(1.D-21), amob0, Ti1
         write (6,*) 'E/n=', E/(1.D-21),'Td'
         write (6,*) 'mui*n=', amob0, '/(m*V*s)'
         write (6,*) 'Ti=', Ti1, 'K'
         print *, ' '
         IF (step.lt.TINY(step)) THEN
            Stop
         END IF
         IF (escale.eq.'log') THEN
            E=dlog(E)+step
            else
            E=E+step
         END IF
         m=m+1
      END DO
      IF (fname.eq.'n') THEN
         PRINT *, ' '
         PRINT *, 'Output saved in file ', filename
         PRINT *, ' '
      END IF
      END ! PROGRAM
C*************************************************************
C QIA1 is the energy-dependent cross section for momentum transfer
C in collisions of ions with neutral particles (note that
C QIA1 is defined here for the particular case of charge exchange)
C AA is a in Eq. (31) of the above mentioned paper
C BB is b in Eq. (31) of the above mentioned paper
      FUNCTION QIA1(Tia,Y)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON AA,BB,ama,ami
      QIA1=2.*(AA-BB*DLOG(8.6183D-5*(Tia)*Y**2))**2
      RETURN
      END FUNCTION
C**************************************************************
C Function QIA11 is the conventional average cross section for momentum
C transfer evaluated at the temperature Ti
      FUNCTION QIA11(A,B,Ti)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      FUN1(Y)=Y**5*DEXP(-Y**2)
      FUN(Y)=FUN1(Y)*QIA1(Ti,Y)
      H=(B-A)/5000
      Y1=A+H
      SOM=FUN(A)+4*FUN(B-H)+FUN(B)
      N=2
      DO WHILE (N.LT.5000)
		SOM1 = SOM + 4*FUN(Y1) + 2*FUN(Y1+H)
          	Y1 = Y1 + 2*H
                N = N + 2
                SOM=SOM1
      END DO
      QIA11=1.d-20/3.*H*SOM1
      RETURN
      END ! FUNCTION
C**************************************************************
C Weighted cross section QIAM
      FUNCTION QIAM(AM,A,B,Tia)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      F1(Y)=Y**5*DEXP(-Y**2)
      F(Y)=F1(Y)*QIA1(Tia,Y)
      FM(Y,AM)=3*(2*AM*Y*DCOSH(2*AM*Y)-DSINH(2*AM*Y))/(2*AM*Y)**3
      H=(B-A)/5000
      Y1=A+H
      SOM = FM(A,AM)*F(A) + 4*FM(B-H,AM)*F(B-H) + FM(B,AM)*F(B)
      N=2
      DO WHILE (N.LT.5000)
		SOM1 = SOM + 4*FM(Y1,AM)*F(Y1) + 2*FM(Y1+H,AM)*F(Y1+H)
		Y1 = Y1 + 2*H
		N = N + 2
	        SOM=SOM1
      END DO
      QIAM=EXP(-AM**2)*(1.D-20/3)*H*SOM1
      RETURN
      END ! FUNCTION
C****************************************************************
C Weighted cross section QIAE
      FUNCTION QIAE(AM,A,B,Tia)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      F1(Y)=Y**5*DEXP(-Y**2)
      F(Y)=F1(Y)*QIA1(Tia,Y)
      FE(Y,AM)=DSINH(2*AM*Y)/(2*AM*Y)
      H=(B-A)/5000
      Y1=A+H
      SOM = FE(A,AM)*F(A) + 4*FE(B-H,AM)*F(B-H) + FE(B,AM)*F(B)
      N=2
      DO WHILE (N.LT.5000)
		SOM1 = SOM + 4*FE(Y1,AM)*F(Y1) + 2*FE(Y1+H,AM)*F(Y1)
                Y1 = Y1 + 2*H
                N = N + 2
                SOM=SOM1
      END DO
      QIAE=EXP(-AM**2)*(1.D-20/3)*H*SOM1
      RETURN
      END ! FUNCTION
C*****************************************************************
      SUBROUTINE GASIONID(GAS,ions)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON AA,BB,ama,ami
      CHARACTER*2 GAS
      CHARACTER*3 ions

      IF (GAS.EQ.'Ar'.and.ions.eq.'Ar+') THEN
        ama=6.6335259d-26
        ami=6.6334348d-26
        AA=7.0
        BB=0.6
      END IF

      IF (GAS.EQ.'He'.and.ions.eq.'He+') THEN
        ama=6.646478151d-27
        ami=6.645567212d-27
        AA=4.88
        BB=0.299
      END IF

      IF (GAS.EQ.'Ne'.and.ions.eq.'Ne+') THEN
        ama=3.350970097d-26
        ami=3.350879003d-26
        AA=5.7
        BB=0.37
      END IF

      IF (GAS.EQ.'Kr'.and.ions.eq.'Kr+') THEN
        ama=1.391532676d-25
        ami=1.391523567d-25
        AA=8.2
        BB=0.57
      END IF

      IF (GAS.EQ.'Xe'.and.ions.eq.'Xe+') THEN
        ama=2.180289265d-25
        ami=2.180280156d-25
        AA=9.2
        BB=0.6
      END IF

      IF (GAS.EQ.'O2'.and.ions.eq.'O2+') THEN
        ama=5.313529333d-26
        ami=5.313438239d-26
        AA=6.2
        BB=0.63
      END IF

      IF (GAS.EQ.'N2'.and.ions.eq.'N2+') THEN
        ama=4.651737647d-26
        ami=4.651646553d-26
        AA=6.5
        BB=0.45
      END IF

      RETURN
      END

