Commit 6afbb16eea83b7812e398d65a3981c9cbdc43ad5
1 parent
26c292bb66
Exists in
master
and in
3 other branches
Correction d'un bug de condition dans ze1
git-svn-id: https://lxsd.femto-st.fr/svn/fvn@56 b657c933-2333-4658-acf2-d3c7c2708721
Showing 1 changed file with 1 additions and 2 deletions Inline Diff
fvn_fnlib/ze1.f90
| function ze1(z) | 1 | 1 | function ze1(z) | |
| ! | 2 | 2 | ! | |
| ! ==================================================== | 3 | 3 | ! ==================================================== | |
| ! Purpose: Compute complex exponential integral E1(z) | 4 | 4 | ! Purpose: Compute complex exponential integral E1(z) | |
| ! Input : z --- Argument of E1(z) | 5 | 5 | ! Input : z --- Argument of E1(z) | |
| ! Output: CE1 --- E1(z) | 6 | 6 | ! Output: CE1 --- E1(z) | |
| ! ==================================================== | 7 | 7 | ! ==================================================== | |
| ! | 8 | 8 | ! | |
| ! Déclaration des variables en passage de paramètre | 9 | 9 | ! Déclaration des variables en passage de paramètre | |
| ! | 10 | 10 | ! | |
| implicit none | 11 | 11 | implicit none | |
| complex(8), intent(in) :: z | 12 | 12 | complex(8), intent(in) :: z | |
| complex(8) :: ze1 | 13 | 13 | complex(8) :: ze1 | |
| ! | 14 | 14 | ! | |
| ! Déclaration des variables locales | 15 | 15 | ! Déclaration des variables locales | |
| ! | 16 | 16 | ! | |
| integer(4) :: k | 17 | 17 | integer(4) :: k | |
| real(8) :: pi,el,x,a0 | 18 | 18 | real(8) :: pi,el,x,a0 | |
| complex(8) :: cr,ct0,ct | 19 | 19 | complex(8) :: cr,ct0,ct | |
| parameter(pi=3.141592653589793D0,el=0.5772156649015328D0) | 20 | 20 | parameter(pi=3.141592653589793D0,el=0.5772156649015328D0) | |
| ! | 21 | 21 | ! | |
| ! traitement en fonction des différents cas | 22 | 22 | ! traitement en fonction des différents cas | |
| ! - Z nul entraîne E1 infini | 23 | 23 | ! - Z nul entraîne E1 infini | |
| ! - module de Z inférieur à 10 ou 20 : formule log+gam+somme | 24 | 24 | ! - module de Z inférieur à 10 ou 20 : formule log+gam+somme | |
| ! - module de Z supérieur à 10 ou 20 : formule asymptotique | 25 | 25 | ! - module de Z supérieur à 10 ou 20 : formule asymptotique | |
| ! | 26 | 26 | ! | |
| x=dreal(z) | 27 | 27 | x=dreal(z) | |
| a0=cdabs(z) | 28 | 28 | a0=cdabs(z) | |
| if (a0==0.0D0) then | 29 | 29 | if (a0==0.0D0) then | |
| ze1 = dcmplx(1.0D+300,0.0D0) | 30 | 30 | ze1 = dcmplx(1.0D+300,0.0D0) | |
| else if ((a0 <= 10.D0).or.(x <= 0.D0.and.a0 <= 20.D0)) then | 31 | 31 | else if ((a0 <= 10.D0).or.(x <= 0.D0.and.a0 <= 20.D0)) then | |
| ze1 = dcmplx(1.0D0,0.0D0) | 32 | 32 | ze1 = dcmplx(1.0D0,0.0D0) | |
| cr = dcmplx(1.0D0,0.0D0) | 33 | 33 | cr = dcmplx(1.0D0,0.0D0) | |
| k=0 | 34 | 34 | k=0 | |
| do while (cdabs(cr)<=cdabs(ze1)*1.0D-15) | 35 | 35 | do while (cdabs(cr)>=cdabs(ze1)*1.0D-15) | |
| k = k+1 | 36 | 36 | k = k+1 | |
| cr = -cr*k*z/(k+1.0D0)**2 | 37 | 37 | cr = -cr*k*z/(k+1.0D0)**2 | |
| ze1 = ze1+cr | 38 | 38 | ze1 = ze1+cr | |
| end do | 39 | 39 | end do | |
| ze1 = -el-cdlog(z)+z*ze1 | 40 | 40 | ze1 = -el-cdlog(z)+z*ze1 | |
| else | 41 | 41 | else | |
| ct0 = dcmplx(0.0D0,0.0D0) | 42 | 42 | ct0 = dcmplx(0.0D0,0.0D0) | |
| do k=120,1,-1 | 43 | 43 | do k=120,1,-1 | |
| ct0 = k/(1.0D0+k/(z+ct0)) | 44 | 44 | ct0 = k/(1.0D0+k/(z+ct0)) | |
| end do | 45 | 45 | end do | |
| ct = 1.0D0/(z+ct0) | 46 | 46 | ct = 1.0D0/(z+ct0) | |
| ze1 = cdexp(-z)*ct | 47 | 47 | ze1 = cdexp(-z)*ct | |
| if (x <= 0.D0.AND.dimag(z) == 0.0) ze1 = ze1-pi*dcmplx(0.D0,1.D0) | 48 | 48 | if (x <= 0.D0.AND.dimag(z) == 0.0) ze1 = ze1-pi*dcmplx(0.D0,1.D0) | |
| end if | 49 | 49 | end if | |
| ! | 50 | 50 | ! | |
| return | 51 | 51 | return | |
| end function | 52 | 52 | end function | |
| 53 | 53 |