Skip to content

Commit

Permalink
version 1.0-2
Browse files Browse the repository at this point in the history
  • Loading branch information
Guy Nason authored and gaborcsardi committed Mar 17, 2010
1 parent e852969 commit 88ccb3b
Show file tree
Hide file tree
Showing 8 changed files with 192 additions and 47 deletions.
24 changes: 15 additions & 9 deletions DESCRIPTION
@@ -1,12 +1,18 @@
Package: NORMT3
Title: Evaluates complex erf, erfc and density of sum of Gaussian and Student's t
Version: 1.0-1
Date: 209-02-137
Title: Evaluates complex erf, erfc, Faddeeva, and density of sum of
Gaussian and Student's t
Version: 1.0-2
Date: 17/03/2010
Author: Guy Nason <g.p.nason@bristol.ac.uk>
Maintainer: Guy Nason <g.p.nason@bristol.ac.uk>
Description: Evaluates the probability density function of the sum of the
Gaussian and Student's t density on 3 degrees of freedom.
Evaluates the p.d.f. of the sphered Student's t density function.
Also evaluates the erf and erfc functions on complex-valued arguments.
License: GPL (>= 2)
Packaged: Fri Feb 13 14:50:55 2009; ripley
Depends: R (>= 2.0)
Description: Evaluates the probability density function of the sum of
the Gaussian and Student's t density on 3 degrees of freedom.
Evaluates the p.d.f. of the sphered Student's t density
function. Also evaluates the erf, and erfc functions on
complex-valued arguments. Thanks to Krishna Myneni the function
is calculates the Faddeeva function also!
License: GPL-2
Packaged: 2012-10-29 08:57:20 UTC; ripley
Repository: CRAN
Date/Publication: 2012-10-29 08:57:20
22 changes: 22 additions & 0 deletions MD5
@@ -0,0 +1,22 @@
fd5fe172ef3ae16f4cc83501ea3215dd *DESCRIPTION
9fdfeb11f7ddba6f3aab9aaa6b1ac5dd *NAMESPACE
adef3a419005deb8c91bade7964cb31c *R/dnormt3.R
dce3e8bb5f353fda3d700bdba01f4651 *R/dst.R
33717d45906e412719aa2035088377f4 *R/erf.R
2b2c7491f277ea56529f7dc734267509 *R/erfc.R
b9b9b1afbe458bdac0b8e34ad1e7163e *R/firstlib.R
411d2dcce18a3bef73c4577aea757ed2 *R/ic1.R
a96ad4f75bd8fb4c9386afb8366caa78 *R/is1.R
2544c29d1fdb284e4b8e0b7faa257080 *R/normt3ip.R
ea867ba161640aaa7004b538e0982f03 *R/wofz.R
69ea61dec509d31c2a8954e2e56596a6 *man/dnormt3.Rd
e607dd933db78b0e85db55bf365174f5 *man/dst.Rd
90037f89dc17dab9a589590bc90a910c *man/erf.Rd
48cb03b2e588846b870566ec6b8ed0f7 *man/erfc.Rd
f28360e4fcdb9b5801cd25324b6dd2b1 *man/ic1.Rd
335b5adaf6aff4553056ca6fbf699dca *man/is1.Rd
b6f6de285b52f88e67c8f850033e3768 *man/normt3ip.Rd
11b45690e9d2298505aba5446c28f75f *man/wofz.Rd
46c60069216c8375f3caa9b8f7ef8b6e *src/IPerfcvec.c
3897cb5687a7258aa7f64e824657a414 *src/IPwofzvec.c
6b033d6d3f308a67c31739db00f4d6de *src/toms680-1.f
5 changes: 5 additions & 0 deletions NAMESPACE
@@ -0,0 +1,5 @@
# Default NAMESPACE created by R
# Remove the previous line if you edit this file

# Export all names
exportPattern(".")
17 changes: 17 additions & 0 deletions R/wofz.R
@@ -0,0 +1,17 @@
"wofz" <-
function(z){

ans <- .C("IPwofzvec",
x=as.double(Re(z)),
y=as.double(Im(z)),
ansx=as.double(Re(z)),
ansy=as.double(Im(z)),
n = as.integer(length(z)),
error = as.integer(0),
PACKAGE="NORMT3")
if (ans$error != 0)
stop(paste("Error code from TOMS 680 was ", ans$error))

return(complex(real=ans$ansx, im=ans$ansy))

}
10 changes: 0 additions & 10 deletions README

This file was deleted.

48 changes: 48 additions & 0 deletions man/wofz.Rd
@@ -0,0 +1,48 @@
\name{wofz}
\alias{wofz}
\title{Faddeeva function}
\description{
Computes the Faddeeva function of a complex valued argument. This function is
\deqn{\exp{-z^2}*erfc(-iz)} .
}
\usage{
wofz(z)
}
\arguments{
\item{z}{Argument of Faddeeva function}
}
\details{
Computes the Faddeeva function of a complex valued argument. This function is
\deqn{\exp{-z^2}*erfc(-iz)}

This function calls FORTRAN code (algorithm TOMS 680)
which computes the Faddeeva function.
}
\value{
The Faddeeva function, w(z)
}
\references{
Poppe, G.P.M. and Wijers, C.M.J. (1990) More efficient computation of the
complex error function. \emph{ACM Transactions on Mathematical
Software}, \bold{16}, 38--46.
}
\author{Krishna Myneni, krishna.myneni@ccreweb.org}


\seealso{\code{\link{erf} \link{erfc}}}
\examples{
options(digits=15)
wofz(0)
#
# Should give 1+0i
#
wofz(1)
#
# Should give 0.367879441171442+0.607157705841394i
#
wofz(complex(re=1, im=1))
#
# Should give 0.304744205256913+0.208218938202832i
#
}
\keyword{math}
54 changes: 54 additions & 0 deletions src/IPwofzvec.c
@@ -0,0 +1,54 @@
#include <math.h>

/*
* IPwofzvec - compute Faddeeva function, w(z) on a vector of complex values
*
*
* Input: x + iy (x and y are double vectors containing the real
and imaginary parts of the complex vector that
you wish to transform)
* Output: ansx + i ansy (ansx and ansy are double vectors containing
the real and imaginary parts of the transformed vector)
* Other variables:
*
* n: Length of each of the vectors x, y, ansx, ansy
*
* error: Error code from the routine that computes Faddeva's function
(computed by TOMS routine 680). Zero means no error.
*/

void IPwofzvec(double *x, double *y, double *ansx, double *ansy, int *n,
int *error)
{
int i, flag;
void IPwofz();
void wofz_();

for(i=0; i< *n; ++i) {
wofz_( x+i, y+i, ansx+i, ansy+i, &flag);

if (flag != 0) {
*error = flag;
return;
}
}
}

void IPwofz(double *a, double *b, double *outa, double *outb, int *error)
{

double d, e;
void wofz_();

*error =0l;

/* Want to work out w(ia+b) */

wofz_( a, b, &d, &e, error);
if (*error) return;

*outa = d;
*outb = e;
}
59 changes: 31 additions & 28 deletions src/toms680.f → src/toms680-1.f
Expand Up @@ -49,6 +49,9 @@ SUBROUTINE WOFZ (XI, YI, U, V, FLAG)
C REFERENCE - GPM POPPE, CMJ WIJERS; MORE EFFICIENT COMPUTATION OF
C THE COMPLEX ERROR-FUNCTION, ACM TRANS. MATH. SOFTWARE.
C
C Revised: 20 August 2009, K. Myneni, promote all literal constants to
C double precision, to avoid loss of accuracy with gfortran.
C
*
*
*
Expand All @@ -65,19 +68,19 @@ SUBROUTINE WOFZ (XI, YI, U, V, FLAG)
*
XABS = DABS(XI)
YABS = DABS(YI)
X = XABS/6.3
Y = YABS/4.4
X = XABS/6.3D0
Y = YABS/4.4D0
*
C
C THE FOLLOWING IF-STATEMENT PROTECTS
C QRHO = (X**2 + Y**2) AGAINST OVERFLOW
C
IF ((XABS.GT.RMAXREAL).OR.(YABS.GT.RMAXREAL)) GOTO 100
*
QRHO = X**2 + Y**2
QRHO = X*X + Y*Y
*
XABSQ = XABS**2
XQUAD = XABSQ - YABS**2
XABSQ = XABS*XABS
XQUAD = XABSQ - YABS*YABS
YQUAD = 2*XABS*YABS
*
A = QRHO.LT.0.085264D0
Expand All @@ -89,18 +92,18 @@ SUBROUTINE WOFZ (XI, YI, U, V, FLAG)
C N IS THE MINIMUM NUMBER OF TERMS NEEDED TO OBTAIN THE REQUIRED
C ACCURACY
C
QRHO = (1-0.85*Y)*DSQRT(QRHO)
N = IDNINT(6 + 72*QRHO)
QRHO = (1.0D0 - 0.85D0*Y)*DSQRT(QRHO)
N = IDNINT(6.0D0 + 72.0D0*QRHO)
J = 2*N+1
XSUM = 1.0/J
XSUM = 1.0D0/J
YSUM = 0.0D0
DO 10 I=N, 1, -1
J = J - 2
XAUX = (XSUM*XQUAD - YSUM*YQUAD)/I
YSUM = (XSUM*YQUAD + YSUM*XQUAD)/I
XSUM = XAUX + 1.0/J
XSUM = XAUX + 1.0D0/J
10 CONTINUE
U1 = -FACTOR*(XSUM*YABS + YSUM*XABS) + 1.0
U1 = -FACTOR*(XSUM*YABS + YSUM*XABS) + 1.0D0
V1 = FACTOR*(XSUM*XABS - YSUM*YABS)
DAUX = DEXP(-XQUAD)
U2 = DAUX*DCOS(YQUAD)
Expand Down Expand Up @@ -129,29 +132,29 @@ SUBROUTINE WOFZ (XI, YI, U, V, FLAG)
H = 0.0D0
KAPN = 0
QRHO = DSQRT(QRHO)
NU = IDINT(3 + (1442/(26*QRHO+77)))
NU = IDINT(3.0D0 + (1442.0D0/(26.0D0*QRHO+77.0D0)))
ELSE
QRHO = (1-Y)*DSQRT(1-QRHO)
H = 1.88*QRHO
H2 = 2*H
KAPN = IDNINT(7 + 34*QRHO)
NU = IDNINT(16 + 26*QRHO)
H = 1.88D0*QRHO
H2 = 2.0D0*H
KAPN = IDNINT(7.0D0 + 34.0D0*QRHO)
NU = IDNINT(16.0D0 + 26.0D0*QRHO)
ENDIF
*
B = (H.GT.0.0)
*
IF (B) QLAMBDA = H2**KAPN
*
RX = 0.0
RY = 0.0
SX = 0.0
SY = 0.0
RX = 0.0D0
RY = 0.0D0
SX = 0.0D0
SY = 0.0D0
*
DO 11 N=NU, 0, -1
NP1 = N + 1
TX = YABS + H + NP1*RX
TY = XABS - NP1*RY
C = 0.5/(TX**2 + TY**2)
C = 0.5D0/(TX*TX + TY*TY)
RX = C*TX
RY = C*TY
IF ((B).AND.(N.LE.KAPN)) THEN
Expand All @@ -162,15 +165,15 @@ SUBROUTINE WOFZ (XI, YI, U, V, FLAG)
ENDIF
11 CONTINUE
*
IF (H.EQ.0.0) THEN
IF (H .EQ. 0.0D0) THEN
U = FACTOR*RX
V = FACTOR*RY
ELSE
U = FACTOR*SX
V = FACTOR*SY
END IF
*
IF (YABS.EQ.0.0) U = DEXP(-XABS**2)
IF (YABS .EQ. 0.0D0) U = DEXP(-XABS*XABS)
*
END IF
*
Expand All @@ -179,11 +182,11 @@ SUBROUTINE WOFZ (XI, YI, U, V, FLAG)
C EVALUATION OF W(Z) IN THE OTHER QUADRANTS
C
*
IF (YI.LT.0.0) THEN
IF (YI .LT. 0.0D0) THEN
*
IF (A) THEN
U2 = 2*U2
V2 = 2*V2
U2 = 2.0D0*U2
V2 = 2.0D0*V2
ELSE
XQUAD = -XQUAD
*
Expand All @@ -194,16 +197,16 @@ SUBROUTINE WOFZ (XI, YI, U, V, FLAG)
IF ((YQUAD.GT.RMAXGONI).OR.
* (XQUAD.GT.RMAXEXP)) GOTO 100
*
W1 = 2*DEXP(XQUAD)
W1 = 2.0D0*DEXP(XQUAD)
U2 = W1*DCOS(YQUAD)
V2 = -W1*DSIN(YQUAD)
END IF
*
U = U2 - U
V = V2 - V
IF (XI.GT.0.0) V = -V
IF (XI .GT. 0.0D0) V = -V
ELSE
IF (XI.LT.0.0) V = -V
IF (XI .LT. 0.0D0) V = -V
END IF
*
RETURN
Expand Down

0 comments on commit 88ccb3b

Please sign in to comment.