By using this site, you agree to our updated Privacy Policy and our Terms of Use. Manage your Cookies Settings.
440,581 Members | 1,966 Online
Bytes IT Community
+ Ask a Question
Need help? Post your question and get tips & solutions from a community of 440,581 IT Pros & Developers. It's quick & easy.

implementation of the complex airy function

P: n/a
I am looking for an implementation of the complex airy function
in C/C++/C#. Can anybody help me find a piece of code for this
function ?

Thanks in advance

Laurent
Jul 22 '05 #1
Share this Question
Share on Google+
4 Replies


P: n/a
In <ee**************************@posting.google.com > la************@orange.fr (Laurent) writes:
I am looking for an implementation of the complex airy function
in C/C++/C#. Can anybody help me find a piece of code for this
function ?


What was your C *and* C++ language question?

Dan
--
Dan Pop
DESY Zeuthen, RZ group
Email: Da*****@ifh.de
Jul 22 '05 #2

P: n/a
la************@orange.fr (Laurent) wrote in message news:<ee**************************@posting.google. com>...
I am looking for an implementation of the complex airy function
in C/C++/C#. Can anybody help me find a piece of code for this
function ?


It's kind of off topic in c.l.c++. But the Airy function can
be implemented in terms of a Bessel function. Bessel fcns
should be in any numerical package.
Socks
Jul 22 '05 #3

P: n/a
Salut Laurent!

/* airy.c CCMATH mathematics library source code.
*
* Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
* This code may be redistributed under the terms of the GNU library
* public license (LGPL). ( See the lgpl.license file for details.)
* ------------------------------------------------------------------------
*/
#include <math.h>
double airy(double x,int df)
{
double f,y,a,b,s;
int p;
double u=.258819403792807,v=.355028053887817;
if(x<=1.7 && x>= -6.9){
y=x*x*x/9.;
if(df){
b= -(a=2./3.);
v*=x*x/2.;
u= -u;
}
else{
a= -(b=1./3.);
u*= -x;
}
for(p=1,f=u+v;;++p){
v*=y/(p*(a+=1.));
u*=y/(p*(b+=1.));
f+=(s=u+v);
if(fabs(s)<1.e-14) break;
}
}
else{
s=1./sqrt(v=3.14159265358979);
y=fabs(x);
if(df) s*=pow(y,.25);
else s/=pow(y,.25);
y*=2.*sqrt(y)/3.;
if(x>0.){
a=12./pow(y,.333);
p=a*a;
if(df) a= -7./36.;
else a=5./36.;
b=2.*(p+y);
f=1.;
u=x=0.;
s*=exp(-y)/2.;
for(; p>0 ;--p,b-=2.){
y=(b*f-(p+1)*x)/(p-1+a/p);
x=f;
u+=(f=y);
}
if(df) f*= -s/u;
else f*=s/u;
}
else{
x=y-v/4.;
y*=2.;
b=.5;
f=s;
v=0.;
if(df) a=2./3.;
else a=1./3.;
for(p=1; (u=fabs(s))>1.e-14 ;++p,b+=1.){
s*=(a+b)*(a-b)/(p*y);
if(fabs(s)>=u) break;
if(!(p&1)){
s= -s;
f+=s;
}
else v+=s;
}
if(df) f=f*sin(x)+v*cos(x);
else f=f*cos(x)-v*sin(x);
}
}
return f;
}
----------------------------------------------------------------------------
------------------
CCMATH: A Mathematics Library
-----------------------------
Version 2.2.0

developed and maintained by
Daniel A. Atkinson
<Da****@aol.com>
Jul 22 '05 #4

P: n/a
Hi Jacob

Thanks a lot for replying with a piece of code. Unfortunatly what I am looking
for is this function with a complex number in argument and not a real number.

Feel free if you have any ideas.

Laurent

"jacob navia" <ja***@jacob.remcomp.fr> wrote in message news:<cb**********@news-reader1.wanadoo.fr>...
Salut Laurent!

/* airy.c CCMATH mathematics library source code.
*
* Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
* This code may be redistributed under the terms of the GNU library
* public license (LGPL). ( See the lgpl.license file for details.)
* ------------------------------------------------------------------------
*/
#include <math.h>
double airy(double x,int df)
{
double f,y,a,b,s;
int p;
double u=.258819403792807,v=.355028053887817;
if(x<=1.7 && x>= -6.9){
y=x*x*x/9.;
if(df){
b= -(a=2./3.);
v*=x*x/2.;
u= -u;
}
else{
a= -(b=1./3.);
u*= -x;
}
for(p=1,f=u+v;;++p){
v*=y/(p*(a+=1.));
u*=y/(p*(b+=1.));
f+=(s=u+v);
if(fabs(s)<1.e-14) break;
}
}
else{
s=1./sqrt(v=3.14159265358979);
y=fabs(x);
if(df) s*=pow(y,.25);
else s/=pow(y,.25);
y*=2.*sqrt(y)/3.;
if(x>0.){
a=12./pow(y,.333);
p=a*a;
if(df) a= -7./36.;
else a=5./36.;
b=2.*(p+y);
f=1.;
u=x=0.;
s*=exp(-y)/2.;
for(; p>0 ;--p,b-=2.){
y=(b*f-(p+1)*x)/(p-1+a/p);
x=f;
u+=(f=y);
}
if(df) f*= -s/u;
else f*=s/u;
}
else{
x=y-v/4.;
y*=2.;
b=.5;
f=s;
v=0.;
if(df) a=2./3.;
else a=1./3.;
for(p=1; (u=fabs(s))>1.e-14 ;++p,b+=1.){
s*=(a+b)*(a-b)/(p*y);
if(fabs(s)>=u) break;
if(!(p&1)){
s= -s;
f+=s;
}
else v+=s;
}
if(df) f=f*sin(x)+v*cos(x);
else f=f*cos(x)-v*sin(x);
}
}
return f;
}
----------------------------------------------------------------------------
------------------
CCMATH: A Mathematics Library
-----------------------------
Version 2.2.0

developed and maintained by
Daniel A. Atkinson
<Da****@aol.com>

Jul 22 '05 #5

This discussion thread is closed

Replies have been disabled for this discussion.