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");
}
}
|
