dtorh.f90

Go to the documentation of this file.

!------------------------------------------------------------------------------!
!> Calculation of toroidal functions \f$P_{n-1/2}^m \left(z\right)\f$ and 
!! \f$Q_{n-1/2}^m \left(z\right)\f$.
!------------------------------------------------------------------------------!
!> \note Copied and adapted from the DTORH1 routine by Segura and Gil 
!! \cite Segura2000 .
!------------------------------------------------------------------------------!
module dtorh
#include <PB3D_macros.h>
    use num_vars, only: dp, pi, max_str_ln
    use messages
    use str_utilities
    
    implicit none
    private
    public dtorh1
    
   character(len=max_str_ln) :: err_msg                                         !< error message
 
contains
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !  INPUT :                                                         !
    !    Z        ARGUMENT OF THE FUNCTIONS                            !
    !    M        ORDER OF THE FUNCTIONS                               !
    !    NMAX     MAXIMUM DEGREE OF THE FUNCTIONS :                    !
    !             WE GET  FUNCTIONS OF ALL THE ORDERS BELOW            !
    !             MIN(NEWN,NMAX). NEWN IS DEFINED BELOW .              !
    !  OUTPUT :                                                        !
    !   *IF MODE IS EQUAL TO 0:                                        !
    !    PL(N)                                                         !
    !             THESE VALUES ARE KEPT IN AN ARRAY                    !
    !    QL(N)                                                         !
    !             THESE VALUES ARE KEPT IN AN ARRAY                    !
    !                                                                  !
    !    NEWN     MAXIMUM ORDER OF FUNCTIONS CALCULATED WHEN           !
    !             PL (NMAX+1) IS LARGER THAN 1/TINY                    !
    !             (OVERFLOW LIMIT = 1/TINY, TINY IS DEFINED BELOW)     !
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !> Calculates toroidal harmonics of a fixed order \c m and degrees
    !! up to \c nmax.
    !!
    !! Optionally, the \c mode can be specified to be different from 0 [default]
    !!  - if \c mode=1:
    !!      - The set of functions evaluated is:
    !!        \f[\frac{P_{n-1/2}^m \left(z\right)}{\Gamma(m+1/2)} \ ,\quad 
    !!          \frac{Q_{n-1/2}^m \left(z\right)}{\Gamma(m+1/2)},\f]
    !!        which are respectively stored in the arrays \c pl(n), \c ql(n).
    !!      - newn refers to this new set of functions.
    !!      - Note 1. and note 2. also apply in this case.
    !!  - if \c mode=2:
    !!      - The code performs as for mode 1, but the restriction \f$ z<20 \f$.
    !!      - For the evaluation of the continued fraction is not considered.
    !!
    !! Also, the parameter \c ipre can be used to select a different precision:
    !!  - For \c ipre=1, the precision is \f$10^{-12}\f$, taking
    !!    \f$\epsilon<10^{-12}\f$
    !!  - For \c ipre=2, the precision is \f$10^{-8}\f$, taking
    !!    \f$\epsilon<10^{-8}\f$
    !!
    !! where \f$\epsilon\f$ controls the accuracy.
    !!
    !! \warning Use \c mode 2 only if high \c m 's for \f$z>20\f$ are required.
    !! The evaluation of the cf may fail to converge for too high \c z 's.
    !!
    !! \note
    !!    -# For a precision of \f$10^{-12}\f$, if \f$z>5\f$ and \f$z/m>0.22\f$,
    !!    the code uses a series expansion for \c pl(0).\n
    !!    When \f$z<20\f$ and \f$z/m<0.22\f$ a continued fraction is applied.
    !!    -# For a precision of \f$10^{-8}\f$, if \f$z>5\f$ and \f$z/m>0.12\f$,
    !!    the code uses a series expansion for \c pl(0).\n
    !!    When \f$z<20\f$ and \f$z/m<0.12\f$ a continued fraction is applied.
    !! 
    !! \return ierr
    integer function dtorh1(z,m,nmax,pl,ql,newn,mode,ipre) result(ierr)
        character(*), parameter :: rout_name = 'DTORH1'
        
        ! input / output
        real(dp), intent(in) :: z                                               !< (real) point at which toroidal harmonics are evaluated
        integer, intent(in) :: m                                                !< order of toroidal harmonics (\c m>0)
        integer, intent(in) :: nmax                                             !< maximum degree of toroidal harmonics (\c nmax>0)
        real(dp), intent(inout) :: pl(0:nmax)                                   !< toroidal harmonics of first kind \f$P_{n-1/2}^m\left(z\right)\f$
        real(dp), intent(inout) :: ql(0:nmax)                                   !< toroidal harmonics of second kind \f$Q_{n-1/2}^m\left(z\right)\f$
        integer, intent(inout) :: newn                                          !< maximum order of functions calculated when <tt> pl(nmax+1)>1/tiny </tt> with <tt>tiny=</tt>\f$10^{-290}\f$
        integer, intent(in), optional :: mode                                   !< mode that controls output
        integer, intent(in), optional :: ipre                                   !< precision
        
        ! local variables
        INTEGER nmaxp,mp,np,n,ical,i
        INTEGER :: mode_loc, ipre_loc
        REAL(dp) over,tinysq,qz,pisq,fl,cc,ar,gamma,fc,qdc1,qargu,&
            &ARGU1,DFACQS,FCP,DD,D1,QM0,DFAC3,DFAC4,PL0
        real(dp), allocatable :: pl_loc(:), ql_loc(:)                           ! local copy of pl and ql, possibly larger
        !REAL(dp) DPPI                                                           ! Is never used?
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        ! THE DIMENSION OF QLMM (INTERNAL ARRAY) MUST BE GREATER THAN M    !
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        REAL(dp)  qlmm(0:1001),pr(2)
        real(dp), PARAMETER :: eps=1.e-14_dp, tiny=1.e-290_dp
        
        ! initialize ierr
        ierr = 0
        
        ! test
        if (nmax.lt.0) then
            ierr = 1
            err_msg = 'Lowest possible degree is 0. Use recursive properties &
                &of toroidal functions if you need lower.'
            chckerr(err_msg)
        end if
        
        ! set optional parameters
        mode_loc = 0
        if (present(mode)) mode_loc = mode
        ipre_loc = 1
        if (present(ipre)) ipre_loc = ipre
        
        over=1._dp/tiny
        tinysq=dsqrt(tiny)
        IF ((ipre_loc.NE.1).AND.(ipre_loc.NE.2)) THEN
            ierr = 1
            err_msg = 'IPRE MUST BE 1 OR 2'
            chckerr(err_msg)
        END IF
        
        pr(1)=.22_dp
        pr(2)=.12_dp
        nmaxp=nmax-1                                                            ! decremented by one to get more logical ouptut
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        !   EPS: REQUIRED ACCURACY FOR THE CONTINUED FRACTION     !
        !        (MODIFIED LENTZ)                                 !
        !   TINY: SMALL PARAMETER TO PREVENT OVERFLOWS IN THE CF  !
        !          (CLOSE TO THE UNDERFLOW LIMIT)                 !
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        IF (z.LE.1._dp) THEN
            ierr = 1
            err_msg = 'IMPROPER ARGUMENT. Z MUST BE GREATER THAN 1'
            chckerr(err_msg)
        END IF
        qz=z
        pisq=dsqrt(pi)
        !DPPI=DSQRT(2._dp)/PISQ                                                  ! Is never used?
        fl=m/2._dp
        cc=abs(float(int(fl))-fl)
        IF (cc.LT.0.4_dp) THEN
            ar=1.
        ELSE
            ar=-1.
        END IF
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        !   WE CHOOSE EXPANSION OR CF FOR PL(0) DEPENDING ON THE VALUES !
        !   OF Z,M AND MODE                                             !
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        ical=1
        IF (m.NE.0) THEN
            IF ((z/m).GT.pr(ipre_loc)) ical=2
            IF (z.LT.5._dp) THEN
                ical=1
            ELSE IF (z.GT.20.d0) THEN
                IF ((mode_loc.NE.2).AND.(ical.EQ.1)) ical=0
            END IF
            IF (ical.EQ.0) THEN
                ierr =1
                chckerr('YOU MUST CHOOSE MODE=2')
            END IF
        ELSE
            IF (z.LT.5.d0) THEN
              ical=1
            ELSE
              ical=2
            END IF
        END IF
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        !   WE USE THE CODE IF NMAX IS GREATER THAN OR EQUAL TO 2  !
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        IF(nmaxp.LT.2) nmaxp=2
        
        ! set up local PL, QL
        allocate(pl_loc(0:nmaxp+1))
        allocate(ql_loc(0:nmaxp+1))
        
        IF (mode_loc.EQ.0) THEN
            gamma=gammah(m,over)*ar*pisq  
            IF (abs(gamma).LT.tiny) THEN
                ierr =1
                err_msg = 'M IS TOO LARGE FOR MODE=0, BETTER TRY MODE=1'
                chckerr(err_msg)
            END IF              
        ELSE 
            gamma=ar
        END IF
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        !   WE EVALUATE THE CONTINUED FRACTION  USING     !
        !   LENTZ-THOMPSON                                !
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        ierr = frac(z,m,0,eps,tinysq,fc)
        chckerr('')
        qdc1=qz*qz-1._dp
        qargu=qz/dsqrt(qdc1)
        !DFAC1=DPPI*GAMMA/PI                                                     ! Is never used?
        !DFAC2=GAMMA/DPPI                                                        ! Is never used?
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        !   WE EVALUATE Q_{-1/2},Q^{1}_{-1/2}               !
        !   USING SLATEC ROUTINES FOR ELLIPTIC FUNCTIONS    !
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        argu1=dsqrt(2._dp/(z+1._dp))
        qlmm(0)=argu1*ellip1(argu1)
        qlmm(1)=-1._dp/dsqrt(2._dp*(qz-1._dp))*ellip2(argu1)
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        !   WE APPLY FORWARD RECURRENCE IN M FOR Q'S  !
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        mp=1
        IF (mode_loc.EQ.0) THEN
          1 IF ((mp.LE.m).AND.(abs(qlmm(mp)).LT.over)) THEN
                qlmm(mp+1)=-2._dp*mp*qargu*qlmm(mp)-&
                    &(mp-0.5_dp)*(mp-0.5_dp)*qlmm(mp-1) 
                mp=mp+1
                GOTO 1 
            ENDIF       
            IF ((mp-1).NE.m) THEN
                ierr = 1
                err_msg = 'M IS TOO LARGE FOR MODE=0, BETTER TRY MODE=1'
                chckerr(err_msg)
            END IF              
        ELSE
            qlmm(0)=qlmm(0)/pisq
            qlmm(1)=qlmm(1)*2._dp/pisq
          2 IF ((mp.LE.m).AND.(abs(qlmm(mp)).LT.over)) THEN
                d1=mp+0.5_dp
                qlmm(mp+1)=-2._dp*mp*qargu*qlmm(mp)/d1&
                    &-(mp-0.5_dp)*qlmm(mp-1)/d1               
                mp=mp+1
                GOTO 2
            ENDIF
            IF ((mp-1).NE.m) THEN
                ierr = 1
                err_msg = 'M IS TOO LARGE FOR MODE=1,2'
                chckerr(err_msg)
            END IF              
        END IF
        !NMMAX=M                                                                 ! Is never used?
        dfacqs=-(gamma/qlmm(m+1))*gamma/pi
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        !   EVALUATION OF PL(0)                             !
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        IF (ical.EQ.1) THEN
            !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
            !   WE CALCULATE THE CF FOR P'S                     !
            !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
            ierr = fracps(qargu,m+1,0,eps,tinysq,fcp) 
            chckerr('')
            dd=m+0.5_dp
            IF (mode_loc.NE.0) THEN
                fcp=fcp/dd
                dfacqs=dfacqs/dd
            END IF
            pl_loc(0)=dfacqs/dsqrt(qdc1)/(1._dp-fcp*qlmm(m)/qlmm(m+1))
        ELSE
            CALL expan(z,mode_loc,ipre_loc,over,qargu,m,pl0)
            pl_loc(0)=pl0
        END IF
        qm0=qlmm(m)
        dfac3=(gamma/qm0)*gamma/pi/(0.5_dp-m)
        pl_loc(1)=pl_loc(0)*fc+dfac3
        np=1 
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        !   WE USE THE RECURRENCE RELATIONS FOR P'S   !
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      3 IF ((np.LE.nmaxp).AND.(abs((np-m+0.5_dp)*pl_loc(np)).LT.over)) THEN
            pl_loc(np+1)=(2._dp*np*z*pl_loc(np)-(np+m-0.5_dp)*pl_loc(np-1))/&
                &(np-m+0.5_dp)
            np=np+1
            GOTO 3
        ENDIF      
        nmaxp=np-1
        newn=nmaxp+1                                                            ! incremented by 1 to get more logical output
        dfac4=(factco(nmaxp,pl_loc(nmaxp+1),m)*gamma)*gamma/pi/(nmaxp+m+0.5_dp)
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        !   WE EVALUATE THE CONTINUED FRACTION USING LENTZ-THOMPSON  !
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        ierr = frac(z,m,nmaxp,eps,tinysq,fc)         
        chckerr('')
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        !   EVALUATION OF PL(NMAX+1) AND PL(NMAX) USING   !
        !   THE WRONSKIAN W{PL(NMAX),QL(NMAX)},           !
        !   THE KNOWN VALUES OF PL(NMAX+1) AND PL(NMAX)   !
        !   THE VALUE OF H = QL(NMAX+1)/QL(NMAX)          !
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        ql_loc(nmaxp)=dfac4/(1._dp-fc*pl_loc(nmaxp)/pl_loc(nmaxp+1))
        ql_loc(nmaxp+1)=ql_loc(nmaxp)*fc
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        !   WE USE THE BACKWARD RECURRENCE RELATION   !
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        DO 4 i=1,nmaxp
            np=nmaxp+1-i
            n=np-1
            ql_loc(n)=((np+np)*z*ql_loc(np)-(np-m+0.5_dp)*ql_loc(np+1))/&
                &(np+m-0.5_dp)
      4 CONTINUE
        
        ! copy local PL and QL back
        pl = pl_loc(0:nmax)
        ql = ql_loc(0:nmax)
    end function dtorh1
    
    ! Adapted to 100k iterations in stead of 10k.
    integer function frac(Z,M,NMAX,EPS,TINYSQ,FC) result(ierr)
        character(*), parameter :: rout_name = 'FRAC'
        
        INTEGER m,nmax,n,mm          
        REAL(dp) z,eps,tinysq,fc,dz2,dn0,dn1,dn2,dn3,dn4,b,a,c0,d0,delta
        
        ! initialize ierr
        ierr = 0
        
        n=nmax
        mm=0
        dz2=z+z
        dn0=n+m
        dn1=n+1._dp
        dn2=dn0+0.5_dp
        dn3=dn0-0.5_dp
        dn4=1._dp-m-m
        b=2._dp*dn1*z/dn2 
        a=1._dp
        fc=tinysq
        c0=fc
        d0=0._dp
     81 d0=b+a*d0
        IF (dabs(d0).LT.tinysq) d0=tinysq
        c0=b+a/c0
        IF (dabs(c0).LT.tinysq) c0=tinysq
        d0=1._dp/d0
        delta=c0*d0
        fc=fc*delta
        mm=mm+1
        a=-(1.d0+dn4/(dn3+mm))  
        b=dz2*(dn1+mm)/(dn2+mm)
        IF (mm.LT.1000000) THEN
            IF (dabs(delta-1.d0).GT.eps) GOTO 81
        END IF
        IF (mm.EQ.1000000) then
            ierr =1
            err_msg = 'CF CONVERGENCE FAILS FOR Z = '//trim(r2str(z))
            chckerr(err_msg)
        END IF
    end function frac
 
    integer function fracps(QZ,M,N,EPS,TINYSQ,FC) result(ierr)
        character(*), parameter :: rout_name = 'FRACPS'
        
        INTEGER m,n,mm         
        REAL(dp) qz,eps,tinysq,fc,dn2,dz2,dm,b,a,c0,d0,delta
        
        ! initialize ierr
        ierr = 0
        
        mm=0
        dn2=n*n
        dz2=qz+qz
        dm=m-0.5_dp
        b=dz2*m
        a=dn2-dm*dm
        fc=tinysq
        c0=fc
        d0=0._dp
     82 d0=b+a*d0
        IF (dabs(d0).LT.tinysq) d0=tinysq
        c0=b+a/c0
        IF (dabs(c0).LT.tinysq) c0=tinysq
        d0=1._dp/d0
        delta=c0*d0
        fc=fc*delta
        mm=mm+1
        a=dn2-(mm+dm)*(mm+dm)  
        b=dz2*(mm+m)
        IF (mm.LT.1000000) THEN
            IF (abs(delta-1._dp).GT.eps) GOTO 82
        END IF
        IF (mm.EQ.1000000) then
            ierr =1
            err_msg = 'CF CONVERGENCE FAILS FOR QZ = '//trim(r2str(qz))
            chckerr(err_msg)
        END IF
    end function fracps
 
    REAL(dp) function ellip2(ak)
        REAL(dp) ak
        REAL(dp) q
        
        q=(1._dp-ak)*(1._dp+ak)
        ellip2=drf(q)-(ak)**2*drd(q)/3.d0
    END FUNCTION ellip2
    
    REAL(dp) function ellip1(ak)
        REAL(dp) ak
        
        ellip1=drf((1._dp-ak)*(1._dp+ak))
    END FUNCTION ellip1
    
    
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !  WE USE SLATEC ROUTINES IN THE EVALUATION OF ELLIPTIC INTEGRALS   !
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    REAL(dp) function drf (y)
        REAL(dp) epslon, errtol
        REAL(dp) c1, c2, c3, e2, e3, lamda
        REAL(dp) mu, s, xn, xndev
        REAL(dp) xnroot, y, yn, yndev, ynroot, zn, zndev, znroot
        LOGICAL first
        SAVE errtol,c1,c2,c3,first
        DATA first /.true./
        
        !***FIRST EXECUTABLE STATEMENT  DRF
        IF (first) THEN
            errtol = 1.e-8_dp
            c1 = 1.0_dp/24.0_dp
            c2 = 3.0_dp/44.0_dp
            c3 = 1.0_dp/14.0_dp
        ENDIF
        first = .false.
        drf = 0.0_dp
        xn = 0._dp
        yn = y
        zn = 1._dp
        
     30 mu = (xn+yn+zn)/3.0_dp
        xndev = 2.0_dp - (mu+xn)/mu
        yndev = 2.0_dp - (mu+yn)/mu
        zndev = 2.0_dp - (mu+zn)/mu
        epslon = max(abs(xndev),abs(yndev),abs(zndev))
        IF (epslon.LT.errtol) GO TO 40
        xnroot = sqrt(xn)
        ynroot = sqrt(yn)
        znroot = sqrt(zn)
        lamda = xnroot*(ynroot+znroot) + ynroot*znroot
        xn = (xn+lamda)*0.250_dp
        yn = (yn+lamda)*0.250_dp
        zn = (zn+lamda)*0.250_dp
        GO TO 30
        
     40 e2 = xndev*yndev - zndev*zndev
        e3 = xndev*yndev*zndev
        s  = 1.0_dp + (c1*e2-0.10_dp-c2*e3)*e2 + c3*e3
        drf = s/sqrt(mu)
    END FUNCTION drf
 
 
    REAL(dp) function drd (y)
        REAL(dp) epslon, errtol
        REAL(dp) c1, c2, c3, c4, ea, eb, ec, ed, ef, lamda
        REAL(dp) mu, power4, sigma, s1, s2, xn, xndev
        REAL(dp) xnroot, y, yn, yndev, ynroot,  zn, zndev, znroot    
        LOGICAL first
        SAVE errtol, c1, c2, c3, c4, first
        DATA first /.true./
        
        !***FIRST EXECUTABLE STATEMENT  DRD
        IF (first) THEN
            errtol = 1.e-8_dp
            c1 = 3.0_dp/14.0_dp
            c2 = 1.0_dp/6.0_dp
            c3 = 9.0_dp/22.0_dp
            c4 = 3.0_dp/26.0_dp
        ENDIF
        first = .false.
        
        ! CALL ERROR HANDLER IF NECESSARY.
        
        drd = 0.0_dp
        
        xn = 0._dp
        yn = y
        zn = 1._dp
        sigma = 0.0_dp
        power4 = 1.0_dp
        
     30 mu = (xn+yn+3.0_dp*zn)*0.20_dp
        xndev = (mu-xn)/mu
        yndev = (mu-yn)/mu
        zndev = (mu-zn)/mu
        epslon = max(abs(xndev), abs(yndev), abs(zndev))
        IF (epslon.LT.errtol) GO TO 40
        xnroot = sqrt(xn)
        ynroot = sqrt(yn)
        znroot = sqrt(zn)
        lamda = xnroot*(ynroot+znroot) + ynroot*znroot
        sigma = sigma + power4/(znroot*(zn+lamda))
        power4 = power4*0.250_dp
        xn = (xn+lamda)*0.250_dp
        yn = (yn+lamda)*0.250_dp
        zn = (zn+lamda)*0.250_dp
        GO TO 30
        
     40 ea = xndev*yndev
        eb = zndev*zndev
        ec = ea - eb
        ed = ea - 6.0_dp*eb
        ef = ed + ec + ec
        s1 = ed*(-c1+0.250_dp*c3*ed-1.50_dp*c4*zndev*ef)
        s2 = zndev*(c2*ef+zndev*(-c3*ec+zndev*c4*ea))
        drd = 3.0_dp*sigma + power4*(1.0_dp+s1+s2)/(mu*sqrt(mu))
    END FUNCTION drd
 
    REAL(dp) function gammah(n,over)
        INTEGER n,i,j
        REAL(dp) over
        
        i=n
        j=2*i-1
        gammah=1._dp
     85 IF ((j.GE.1).AND.(gammah.LT.over)) THEN
            gammah=gammah*j/2._dp
            i=i-1
            j=2*i-1
            GOTO 85
        ENDIF  
        IF (j.GT.1) THEN 
            gammah=0._dp
        END IF   
    END FUNCTION gammah
 
    REAL(dp) function factco(n,pl,m)
        INTEGER n,m,j  
        REAL(dp) pl,x1,x2 
        
        factco=1._dp/pl
        x1=m+0.5_dp
        x2=-m+0.5_dp
        j=n
    86  IF (j.GE.0) THEN
            factco=factco*(j+x1)/(j+x2)
            j=j-1
            GOTO 86
        ENDIF
    END FUNCTION factco
 
    SUBROUTINE expan(Z,MODE,IPRE,OVER,QARGU,M,PL)
        INTEGER mode,ipre,m,k,i 
        REAL(dp) z,over,qargu,pl,pisq,db,fl,cc,dz,gamma,ar,&
            &DFAC,DF1,A0,Z2I,DA1,DA2,DELTA,SUM,DCC,DCCP,DKK
        REAL(dp) preci(2)
        
        preci(1)=1.e-13_dp
        preci(2)=1.e-9_dp
        pisq=dsqrt(pi)
        db=2._dp*dlog(2._dp)
        fl=m/2.
        cc=abs(float(int(fl))-fl)
        IF (cc.LT.0.4_dp) THEN
            ar=1.
        ELSE
            ar=-1.
        END IF      
        dz=1._dp
        DO 100 i=1,m
            dz=dz/qargu
    100 CONTINUE
        IF (mode.EQ.0) THEN
            gamma=gammah(m,over)*ar*pisq 
        ELSE 
            gamma=ar
        END IF
        dfac=2._dp/pi*dz*gamma/pisq
        df1=dlog(2._dp*z)
        a0=1._dp/dsqrt(2._dp*z)
        z2i=1._dp/(z*z)
        delta=1._dp
        sum=0._dp
        k=0
        da2=factor(m,0)
        da1=db+psi(m,0)
     87 delta=(df1+da1)*da2*a0
        sum=sum+delta
        dcc=0.5_dp+m+2._dp*k
        dccp=dcc+1._dp
        dkk=k+1._dp
        da2=da2*dccp*dcc/(dkk*dkk)
        da1=da1+1._dp/dkk-1._dp/dccp-1._dp/dcc
        k=k+1
        a0=a0*0.25_dp*z2i
        IF (abs(delta/sum).GT.preci(ipre)) GOTO 87
        pl=sum*dfac
    END SUBROUTINE expan
 
    REAL(dp) function factor(m,k)
        INTEGER m,k,n,i,j
        REAL(dp) x1
        
        factor=1._dp
        IF (k.GE.1) THEN
        x1=m+0.5_dp
        n=2*k-1
        i=k
        j=n
     88 IF (j.GE.0) THEN
            factor=factor*(j+x1)/i/i
            j=j-1
            i=i-1
            IF (i.EQ.0) i=1
                GOTO 88
            ENDIF 
        END IF  
    END FUNCTION factor
 
    REAL(dp) function psi(m,k)
        INTEGER m,k,n,j,i
        REAL(dp) factr1,factr2
        psi=0._dp
        factr1=0._dp
        factr2=0._dp
        n=2*k+m
        j=1
        IF (k.GE.1) THEN
         89 IF (j.LE.k) THEN
                factr1=factr1+1._dp/j
                j=j+1
                GOTO 89
            ENDIF
        END IF
        i=1
     90 IF (i.LE.n) THEN
            factr2=factr2+1._dp/(2._dp*i-1._dp)
            i=i+1
            GOTO 90
        ENDIF
        IF (n.EQ.0) factr2=0._dp
        psi=factr1-2._dp*factr2
    END FUNCTION psi
end module