Contribution de Sylvain :

fonction e1 avec arguments complexes ze1 git-svn-id: https://lxsd.femto-st.fr/svn/fvn@55 b657c933-2333-4658-acf2-d3c7c2708721

Contribution de Sylvain :
fonction e1 avec arguments complexes ze1 git-svn-id: https://lxsd.femto-st.fr/svn/fvn@55 b657c933-2333-4658-acf2-d3c7c2708721
wdaniau
1 parent 6164710c3d
Showing 5 changed files with 2347 additions and 2298 deletions Side-by-side Diff
doc/fvn.pdf
doc/fvn.tex
fvn_fnlib/Makefile
fvn_fnlib/fvn_fnlib.f90
fvn_fnlib/ze1.f90
@@ -66,10 +66,11 @@
     Copyright, this License, and the Availability note are retained,
     and a notice that the code was modified is included.
  
-\subsubsection*{Authors}
-As of the day this manuel is written there's only one author of fvn :\newline
-William Daniau\newline
-william.daniau@femto-st.fr\newline
+\subsubsection*{Authors and Contributors}
+\begin{itemize}
+ \item William Daniau (william.daniau@femto-st.fr)
+ \item Sylvain Ballandras (sylvain.ballandras@femto-st.fr)
+\end{itemize}
  
 \section{Naming scheme and convention}
 The naming scheme of the routines is as follow :
  
  
@@ -1228,15 +1229,28 @@
  e1(x)
 \end{verbatim}
 \begin{itemize}
- \item x is a real
+ \item x is a real or complex
 \end{itemize}
-This function evaluates the exponential integral for argument greater than 0 and the Cauchy principal value for argument less than 0. It is define by equation \ref{e1} for $x \neq 0$.
+For a real argument, this function evaluates the exponential integral for argument greater than 0 and the Cauchy principal value for argument less than 0. It is define by equation \ref{e1} for $x \neq 0$.
 \begin{equation}
 \label{e1}
  e1(x)= \int _{x} ^\infty {e^{-t}\over t}dt
 \end{equation}
  
-Specific interfaces : \verb'e1,de1'
+For a complex argument, the notation in equation \ref{e1cplx} is used (Abramowitz and Stegun, p.228 \url{http://www.math.ucla.edu/~cbm/aands/page_228.htm}):
+\begin{equation}
+\label{e1cplx}
+ e1(z) = \int _{z} ^\infty {e^{-t}\over t}dt \textrm{~with~} \left|  arg(z) \right| < \pi
+\end{equation}
+For positive values of real part of $z$, this can be written as in equation \ref{e1cplx_pos} :
+\begin{equation}
+\label{e1cplx_pos}
+ e1(z)= \int _{1} ^\infty {e^{-tz}\over t}dt \textrm{~with~} Re(z) > 0
+\end{equation}
+
+ 
+
+Specific interfaces : \verb'e1,de1,ze1'
  
 \subsubsection{ali}
 \begin{verbatim}
@@ -61,7 +61,7 @@
 zlngam.o zlnrel.o zlog10.o zpsi.o  \
 zsinh.o ztanh.o ztan.o besyn.o \
 besjn.o dbesyn.o dbesjn.o beskn.o \
-besin.o dbeskn.o dbesin.o
+besin.o dbeskn.o dbesin.o ze1.o
  
 lib:$(objects)
  
@@ -282,6 +282,9 @@
       real(8) function de1(x)
             real(8) :: x
       end function de1
+      complex(8) function ze1(x)
+            complex(8) :: x
+      end function
 end interface e1
  
 !!!!!!!!!!!!!!!
+function  ze1(z)
+!
+!       ====================================================
+!       Purpose: Compute complex exponential integral E1(z)
+!       Input :  z   --- Argument of E1(z)
+!       Output:  CE1 --- E1(z)
+!       ====================================================
+!
+! Déclaration des variables en passage de paramètre
+!
+implicit none
+complex(8), intent(in) :: z
+complex(8) :: ze1
+!
+!  Déclaration des variables locales
+!
+integer(4) :: k
+real(8) :: pi,el,x,a0
+complex(8) :: cr,ct0,ct
+parameter(pi=3.141592653589793D0,el=0.5772156649015328D0)
+!
+! traitement en fonction des différents cas
+! - Z nul entraîne E1 infini
+! - module de Z inférieur à 10 ou 20 : formule log+gam+somme
+! - module de Z supérieur à 10 ou 20 : formule asymptotique
+!
+ x=dreal(z)
+ a0=cdabs(z)
+ if (a0==0.0D0) then
+   ze1 = dcmplx(1.0D+300,0.0D0)
+ else if ((a0 <= 10.D0).or.(x <= 0.D0.and.a0 <= 20.D0)) then
+   ze1 = dcmplx(1.0D0,0.0D0)
+   cr = dcmplx(1.0D0,0.0D0)
+   k=0
+   do while (cdabs(cr)<=cdabs(ze1)*1.0D-15)
+     k = k+1
+     cr = -cr*k*z/(k+1.0D0)**2
+     ze1 = ze1+cr
+   end do
+   ze1 = -el-cdlog(z)+z*ze1
+ else
+   ct0 = dcmplx(0.0D0,0.0D0)
+   do k=120,1,-1
+      ct0 = k/(1.0D0+k/(z+ct0))
+   end do
+   ct = 1.0D0/(z+ct0)
+   ze1 = cdexp(-z)*ct
+   if (x <= 0.D0.AND.dimag(z) == 0.0) ze1 = ze1-pi*dcmplx(0.D0,1.D0)
+ end if
+!
+ return
+end function
...	...	@@ -66,10 +66,11 @@
66	66	Copyright, this License, and the Availability note are retained,
67	67	and a notice that the code was modified is included.
68	68
69		-\subsubsection*{Authors}
70		-As of the day this manuel is written there's only one author of fvn :\newline
71		-William Daniau\newline
72		-william.daniau@femto-st.fr\newline
	69	+\subsubsection*{Authors and Contributors}
	70	+\begin{itemize}
	71	+ \item William Daniau (william.daniau@femto-st.fr)
	72	+ \item Sylvain Ballandras (sylvain.ballandras@femto-st.fr)
	73	+\end{itemize}
73	74
74	75	\section{Naming scheme and convention}
75	76	The naming scheme of the routines is as follow :
76	77
77	78
...	...	@@ -1228,15 +1229,28 @@
1228	1229	e1(x)
1229	1230	\end{verbatim}
1230	1231	\begin{itemize}
1231		- \item x is a real
	1232	+ \item x is a real or complex
1232	1233	\end{itemize}
1233		-This function evaluates the exponential integral for argument greater than 0 and the Cauchy principal value for argument less than 0. It is define by equation \ref{e1} for $x \neq 0$.
	1234	+For a real argument, this function evaluates the exponential integral for argument greater than 0 and the Cauchy principal value for argument less than 0. It is define by equation \ref{e1} for $x \neq 0$.
1234	1235	\begin{equation}
1235	1236	\label{e1}
1236	1237	e1(x)= \int _{x} ^\infty {e^{-t}\over t}dt
1237	1238	\end{equation}
1238	1239
1239		-Specific interfaces : \verb'e1,de1'
	1240	+For a complex argument, the notation in equation \ref{e1cplx} is used (Abramowitz and Stegun, p.228 \url{http://www.math.ucla.edu/~cbm/aands/page_228.htm}):
	1241	+\begin{equation}
	1242	+\label{e1cplx}
	1243	+ e1(z) = \int _{z} ^\infty {e^{-t}\over t}dt \textrm{~with~} \left\| arg(z) \right\| < \pi
	1244	+\end{equation}
	1245	+For positive values of real part of $z$, this can be written as in equation \ref{e1cplx_pos} :
	1246	+\begin{equation}
	1247	+\label{e1cplx_pos}
	1248	+ e1(z)= \int _{1} ^\infty {e^{-tz}\over t}dt \textrm{~with~} Re(z) > 0
	1249	+\end{equation}
	1250	+
	1251	+
	1252	+
	1253	+Specific interfaces : \verb'e1,de1,ze1'
1240	1254
1241	1255	\subsubsection{ali}
1242	1256	\begin{verbatim}
...	...	@@ -61,7 +61,7 @@
61	61	zlngam.o zlnrel.o zlog10.o zpsi.o \
62	62	zsinh.o ztanh.o ztan.o besyn.o \
63	63	besjn.o dbesyn.o dbesjn.o beskn.o \
64		-besin.o dbeskn.o dbesin.o
	64	+besin.o dbeskn.o dbesin.o ze1.o
65	65
66	66	lib:$(objects)
67	67
...	...	@@ -282,6 +282,9 @@
282	282	real(8) function de1(x)
283	283	real(8) :: x
284	284	end function de1
	285	+ complex(8) function ze1(x)
	286	+ complex(8) :: x
	287	+ end function
285	288	end interface e1
286	289
287	290	!!!!!!!!!!!!!!!
	1	+function ze1(z)
	2	+!
	3	+! ====================================================
	4	+! Purpose: Compute complex exponential integral E1(z)
	5	+! Input : z --- Argument of E1(z)
	6	+! Output: CE1 --- E1(z)
	7	+! ====================================================
	8	+!
	9	+! Déclaration des variables en passage de paramètre
	10	+!
	11	+implicit none
	12	+complex(8), intent(in) :: z
	13	+complex(8) :: ze1
	14	+!
	15	+! Déclaration des variables locales
	16	+!
	17	+integer(4) :: k
	18	+real(8) :: pi,el,x,a0
	19	+complex(8) :: cr,ct0,ct
	20	+parameter(pi=3.141592653589793D0,el=0.5772156649015328D0)
	21	+!
	22	+! traitement en fonction des différents cas
	23	+! - Z nul entraîne E1 infini
	24	+! - module de Z inférieur à 10 ou 20 : formule log+gam+somme
	25	+! - module de Z supérieur à 10 ou 20 : formule asymptotique
	26	+!
	27	+ x=dreal(z)
	28	+ a0=cdabs(z)
	29	+ if (a0==0.0D0) then
	30	+ ze1 = dcmplx(1.0D+300,0.0D0)
	31	+ else if ((a0 <= 10.D0).or.(x <= 0.D0.and.a0 <= 20.D0)) then
	32	+ ze1 = dcmplx(1.0D0,0.0D0)
	33	+ cr = dcmplx(1.0D0,0.0D0)
	34	+ k=0
	35	+ do while (cdabs(cr)<=cdabs(ze1)*1.0D-15)
	36	+ k = k+1
	37	+ cr = -crkz/(k+1.0D0)**2
	38	+ ze1 = ze1+cr
	39	+ end do
	40	+ ze1 = -el-cdlog(z)+z*ze1
	41	+ else
	42	+ ct0 = dcmplx(0.0D0,0.0D0)
	43	+ do k=120,1,-1
	44	+ ct0 = k/(1.0D0+k/(z+ct0))
	45	+ end do
	46	+ ct = 1.0D0/(z+ct0)
	47	+ ze1 = cdexp(-z)*ct
	48	+ if (x <= 0.D0.AND.dimag(z) == 0.0) ze1 = ze1-pi*dcmplx(0.D0,1.D0)
	49	+ end if
	50	+!
	51	+ return
	52	+end function