Siguiente Anterior

Integrales de Fresnel

Algunas páginas del Curso Interactivo de Física describen situaciones físicas en las que es necesario calcular las integrales S(t) y C(t).

S(t)= 0 t sin( π 2 u 2 )du C(t)= 0 t cos( π 2 u 2 )du

public class Fresnel {
    static final double EPS=6.0e-8;
    static final int MAXIT=100;
    static final double FPMIN=1.0e-30;
    static final double XMIN=1.5;
//Here EPS is the relative error; MAXIT is the maximum number 
of iterations allowed; FPMIN //is a number near the smallest representable floating-point number;
XMIN is the dividing line //between using the series and continued fraction. //valores de las funciones double s; //integral S double c; //integral C //Computes the Fresnel integrals S(x) and C(x) for all real x. void integral(double x){ int n, k; boolean odd; double a,ax,fact,pix2,sign,sum,sumc,sums,term,test; Complejo b,cc,d,h,del,cs; ax=Math.abs(x); if (ax < Math.sqrt(FPMIN)) { //Secial case: avoid failure of convergence //test because of underflow. *s=0.0; c=ax; } else if (ax <= XMIN) { //Evaluate both series simultaneously. sum=sums=0.0; sumc=ax; sign=1.0; fact=Math.PI*ax*ax/2; odd=true; term=ax; n=3; for (k=1;k<=MAXIT;k++) { term *= fact/k; sum += sign*term/n; test=Math.abs(sum)*EPS; if (odd) { sign = -sign; sums=sum; sum=sumc; } else { sumc=sum; sum=sums; } if (term < test) break; odd=!odd; n += 2; } if (k > MAXIT) nrerror("series failed in frenel"); s=sums; c=sumc; } else { // Evaluate continued fraction by modi_ed Lentz's method (x5.2). pix2=Math.PI*ax*ax; b=new Complejo(1.0,-pix2); cc=new Complejo(1.0/FPMIN,0.0); d=h=Complejo.cociente(new Complejo(1.0, 0.0),b); n = -1; for (k=2;k<=MAXIT;k++) { n += 2; a = -n*(n+1); b=Complejo.suma(b,new Complejo(4.0,0.0)); d=Complejo.cociente(new Complejo(1.0, 0.0),
Complejo.suma(Complejo.producto(a,d),b)); //Denominators cannot be zero. cc=Complejo.suma(b,Complejo.cociente(new Complejo(a,0.0),cc)); del=Complejo.producto(cc,d); h=Complejo.producto(h,del); if (Math.abs(del.real-1.0)+Math.abs(del.imag) < EPS) break; } if (k > MAXIT) nrerror("cf failed in frenel"); h=Complejo.producto(new Complejo(ax,-ax),h); cs=Complejo.producto(new Complejo(0.5,0.5),
Complejo.diferencia(new Complejo(1.0, 0.0),
Complejo.producto(new Complejo(Math.cos(0.5*pix2),Math.sin(0.5*pix2)),h))); c=cs.real; s=cs.imag; } if (x < 0.0) { // Use antisymmetry. c = -c; s = -s; } } void nrerror(String texto){ System.out.println(texto); } } public class Complejo{ double real; double imag; public Complejo() { real=0.0; imag=0.0; } public Complejo(double real, double imag){ this.real=real; this.imag=imag; } public static Complejo conjugado(Complejo c){ return new Complejo(c.real, -c.imag); } public static Complejo opuesto(Complejo c){ return new Complejo(-c.real, -c.imag); } public double modulo(){ return Math.sqrt(real*real+imag*imag); } //devuelve el ángulo en grados public double argumento(){ double angulo=Math.atan2(imag, real);
//devuelve el ángulo entre -PI y +PI if(angulo<0) angulo=2*Math.PI+angulo; return angulo*180/Math.PI; } //suma de dos números complejos public static Complejo suma(Complejo c1, Complejo c2){ double x=c1.real+c2.real; double y=c1.imag+c2.imag; return new Complejo(x, y); } public static Complejo diferencia(Complejo c1, Complejo c2){ double x=c1.real-c2.real; double y=c1.imag-c2.imag; return new Complejo(x, y); } //producto de dos números complejos public static Complejo producto(Complejo c1, Complejo c2){ double x=c1.real*c2.real-c1.imag*c2.imag; double y=c1.real*c2.imag+c1.imag*c2.real; return new Complejo(x, y); } //producto de un complejo por un número real public static Complejo producto(Complejo c, double d){ double x=c.real*d; double y=c.imag*d; return new Complejo(x, y); } //producto de un número real por un complejo public static Complejo producto(double d, Complejo c){ double x=c.real*d; double y=c.imag*d; return new Complejo(x, y); } //cociente de dos números complejos public static Complejo cociente(Complejo c1, Complejo c2){
double aux, x, y;
aux=c2.real*c2.real+c2.imag*c2.imag; x=(c1.real*c2.real+c1.imag*c2.imag)/aux; y=(c1.imag*c2.real-c1.real*c2.imag)/aux; return new Complejo(x, y); } //cociente entre un número complejo y un número real public static Complejo cociente(Complejo c, double d){ double x=c.real/d; double y=c.imag/d; return new Complejo(x, y); } //el número e elevado a un número complejo public static Complejo exponencial(Complejo c){ double x=Math.cos(c.imag)*Math.exp(c.real); double y=Math.sin(c.imag)*Math.exp(c.real); return new Complejo(x, y); } //raíz cuadrada de un número positivo o negativo public static Complejo csqrt(double d){ if(d>=0) return new Complejo(Math.sqrt(d), 0); return new Complejo(0, Math.sqrt(-d)); } //función auxiliar para la potencia de un número complejo private static double potencia(double base, int exponente){ double resultado=1.0; for(int i=0; i<exponente; i++){ resultado*=base; } return resultado; } //función auxiliar para la potencia de un número complejo private static double combinatorio(int m, int n){ long num=1; long den=1; for(int i=m; i>m-n; i--){ num*=i; } for(int i=2; i<=n; i++){ den*=i; } return (double)num/den; } //potencia de un número complejo public static Complejo potencia(Complejo c, int exponente){ double x=0.0, y=0.0; int signo; for(int i=0; i<=exponente; i++){ signo=(i%2==0)?+1:-1; //parte real x+=combinatorio(exponente, 2*i)*potencia(c.real, exponente-2*i)
*potencia(c.imag, 2*i)*signo; if(exponente==2*i) break; //parte imaginaria y+=combinatorio(exponente, 2*i+1)*potencia(c.real, exponente-(2*i+1))
*potencia(c.imag, 2*i+1)*signo; } return new Complejo(x, y); } //traduce un número complejo a un string public String toString(){ if(imag>0) return new String((double)Math.round(100*real)/100+" +
"+(double)Math.round(100*imag)/100+"*i"); return new String((double)Math.round(100*real)/100+" - "
+(double)Math.round(-100*imag)/100+"*i"); } }


Press W. H., Teukolsky S. A., Vetterling W. T., Flannery B. P. Numerical Recipes in C, Second edition, Special functions. Fresnel Integrals, Cosine and Sine Integrals  Chapter 6º. Cambridge University Press. Código en C adaptado por el autor al lenguaje Java.

Siguiente Anterior