lambertw.f90

Go to the documentation of this file.

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! Time-stamp: <2015-05-21 13:35:03 pbrowne>
!!!
!!!    Subroutine to implement the lambertw function
!!!    Copyright (C) 2015 Mengbin Zhu 
!!!
!!!    This program is free software: you can redistribute it and/or modify
!!!    it under the terms of the GNU General Public License as published by
!!!    the Free Software Foundation, either version 3 of the License, or
!!!    (at your option) any later version.
!!!
!!!    This program is distributed in the hope that it will be useful,
!!!    but WITHOUT ANY WARRANTY; without even the implied warranty of
!!!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!!!    GNU General Public License for more details.
!!!
!!!    You should have received a copy of the GNU General Public License
!!!    along with this program.  If not, see <http://www.gnu.org/licenses/>.
!!!
!!!    Email: zhumengbin @ gmail.com  
!!!    Mail:  School of Mathematical and Physical Sciences,
!!!           University of Reading,
!!!       Reading, UK
!!!       RG6 6BB
!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

SUBROUTINE lambertw(K,X,W)

!######################################################################
!Lambert W Function Numerical Solver
!In mathematics, the Lambert W function, also called the omega function 
!or product logarithm, is a set of functions, namely the branches of 
!the inverse relation of the function z = f(W) = We^W 
!where e^W is the exponential function and W is any complex number.
!
!University of Reading, 
!Dept. of Meteorology, 
!Mengbin Zhu 14/04/2015
!######################################################################

IMPLICIT NONE

INTEGER, INTENT(IN)  :: K
REAL(KIND=KIND(1.0d0)),    INTENT(IN)  :: X
REAL(KIND=KIND(1.0d0)),    INTENT(OUT) :: W

INTEGER              :: I,J,Q,R,Q2,Q3

REAL(KIND=KIND(1.0d0)),PARAMETER       :: EPS=4.0E-16
REAL(KIND=KIND(1.0d0)),PARAMETER       :: EM1=0.3678794411714423215955237701614608

REAL(KIND=KIND(1.0d0))                 :: PP,EE,TT2,L1,L2

IF(k==0) THEN     !Lambert W-Function 0 Branch
    IF(x<-em1) THEN
        print*,'Lambert W-Function: Bad Argument! STOP!'
        RETURN
    END IF
    IF(x==0.0) THEN
        w=0.0
        RETURN
    END IF
    IF(x<(-em1+1e-4)) THEN
        q=x+em1
        r=sqrt(REAL(q))
        q2=q*q
        q3=q2*q
        w=-1.0 &
        & +2.331643981597124203363536062168*r    &
        & -1.812187885639363490240191647568*q    &
        & +1.936631114492359755363277457668*r*q  &
        & -2.353551201881614516821543561516*q2   &
        & +3.066858901050631912893148922704*r*q2 &
        & -4.175335600258177138854984177460*q3   &
        & +5.858023729874774148815053846119*r*q3 &
        & -8.401032217523977370984161688514*q3*q
        RETURN
    END IF
    IF(x<1.0) THEN
        pp=sqrt(2.0*(2.7182818284590452353602874713526625*x+1.0))
        w=-1.0+pp*(1.0+pp*(-0.333333333333333333333+pp*0.152777777777777777777777))
    ELSE
        w=log(x)
    END IF
    !IF(X>3.0) THEN
    !    W=W-LOG(W)
    !END IF

    tt2=1.0
    i=0
    DO WHILE (abs(tt2)>(eps*(1.0+abs(w))))
        !PRINT*,"Cannot go out of the cycle 0 Branch."
        i=i+1
        ee=exp(w)
        tt2=w*ee-x
        pp=w+1.0
        tt2=tt2/(ee*pp-0.5*(pp+1.0)*tt2/pp)
        w=w-tt2
        !PRINT*,W
        IF(i==20) THEN
            w = w + 1.0e-6
        END IF
        IF(i==40) EXIT
    END DO

ELSE IF (k==-1) THEN  !Lambert W-Function -1 Branch.
    !IF(X<-EM1 .OR. X>0.0) THEN
    !    PRINT*,'Lambert W-Function -1 Branch: Bad Argument! STOP!'
    !    !STOP
    !    RETURN
    !END IF
    !PRINT*,"THE -1 BRANCH"
    IF(x==0.0) THEN
        w=0.0
    ELSE 
        !PRINT*,X
        pp=1.0
        IF(x<(-1e-6)) THEN
            !PRINT*,"IS THIS REAL?"
            pp=-sqrt(2.0*(2.7182818284590452353602874713526625*x+1.0))
            w=-1.0+pp*(1.0+pp*(-0.333333333333333333333+pp*0.152777777777777777777777))
        ELSE
            l1=log(-x)
            l2=log(-l1)
            w=l1-l2+l2/l1
        END IF
        IF(abs(pp)<1e-4) THEN
            RETURN
        END IF
        tt2=1.0
        i=0
        DO WHILE(abs(tt2)>eps*(1.0+abs(w)))
            !PRINT*,"Cannot go out of the cycle"
            i=i+1
            ee=exp(w)
            tt2=w*ee-x
            pp=w+1.0
            tt2=tt2/(ee*pp-0.5*(pp+1.0)*tt2/pp)
            w=w-tt2
            !PRINT*,W
            IF(i==20) THEN
                w = w + 1.0e-6
            END IF
            IF(i==40) EXIT
        END DO
    ENDIF

ELSE
    print*,"Lambert W-Function Branch Should Be -1 OR 0!"
    !STOP
    RETURN
END IF

END SUBROUTINE lambertw