]> Función de Airy
Anterior

Función de Airy, Ai

public class Numero {
    public double valor;
  public Numero(double x) {
      this.valor=x;
  }
}
public class Airy {
static final double EPS=1.0e-10;
static final double FPMIN=1.0e-30;
static final int MAXIT=10000;
static final double XMIN= 2.0;
static final double PI= 3.141592653589793;


public static  double chebev(double a, double b, double[] c, int m, double x){
        double d=0.0,dd=0.0,sv,y,y2;
        int j;
        if ((x-a)*(x-b) > 0.0) nrerror("x not in range in routine chebev");
        y2=2.0*(y=(2.0*x-a-b)/(b-a));
        for (j=m-1;j>=1;j--) {
          sv=d;
          d=y2*d-dd+c[j];
          dd=sv;
        }
        return y*d-dd+0.5*c[0];
  }
  public static void beschb(double x, Numero gam1, Numero gam2, Numero gampl, 
Numero gammi){ final int NUSE1=7; final int NUSE2=8; double xx; final double[] c1 = {-1.142022680371168e0,6.5165112670737e-3,
3.087090173086e-4,-3.4706269649e-6,6.9437664e-9, 3.67795e-11,-1.356e-13}; final double[] c2 = {1.843740587300905e0,-7.68528408447867e-2,
1.2719271366546e-3,-4.9717367042e-6,-3.31261198e-8, 2.423096e-10,-1.702e-13,-1.49e-15}; xx=8.0*x*x-1.0; gam1.valor=chebev(-1.0,1.0,c1,NUSE1,xx); gam2.valor=chebev(-1.0,1.0,c2,NUSE2,xx); gampl.valor= gam2.valor-x*gam1.valor; gammi.valor= gam2.valor+x*gam1.valor; } public static void bessjy(double x, double xnu, Numero rj, Numero ry, Numero rjp,
Numero ryp) { int i,isign,l,nl; double a,b,br,bi,c,cr,ci,d,del,del1,den,di,dlr,dli,dr,e,f,fact,fact2, fact3,ff,gam,gam1=0.0,gam2=0.0,gammi=0.0,gampl=0.0,h,p,pimu,pimu2,q,r,rjl, rjl1,rjmu,rjp1,rjpl,rjtemp,ry1,rymu,rymup,rytemp,sum,sum1, temp,w,x2,xi,xi2,xmu,xmu2; if (x <= 0.0 || xnu < 0.0) nrerror("bad arguments in bessjy"); nl=(x < XMIN ? (int)(xnu+0.5) : IMAX(0,(int)(xnu-x+1.5))); xmu=xnu-nl; xmu2=xmu*xmu; xi=1.0/x; xi2=2.0*xi; w=xi2/PI; isign=1; h=xnu*xi; if (h < FPMIN) h=FPMIN; b=xi2*xnu; d=0.0; c=h; for (i=1;i<=MAXIT;i++) { b += xi2; d=b-d; if (Math.abs(d) < FPMIN) d=FPMIN; c=b-1.0/c; if (Math.abs(c) < FPMIN) c=FPMIN; d=1.0/d; del=c*d; h=del*h; if (d < 0.0) isign = -isign; if (Math.abs(del-1.0) < EPS) break; } if (i > MAXIT) nrerror("x too large in bessjy; try asymptotic expansion"); rjl=isign*FPMIN; rjpl=h*rjl; rjl1=rjl; rjp1=rjpl; fact=xnu*xi; for (l=nl;l>=1;l--) { rjtemp=fact*rjl+rjpl; fact -= xi; rjpl=fact*rjtemp-rjl; rjl=rjtemp; } if (rjl == 0.0) rjl=EPS; f=rjpl/rjl; if (x < XMIN) { x2=0.5*x; pimu=PI*xmu; fact = (Math.abs(pimu) < EPS ? 1.0 : pimu/Math.sin(pimu)); d = -Math.log(x2); e=xmu*d; fact2 = (Math.abs(e) < EPS ? 1.0 : sinh(e)/e); Numero g1=new Numero(gam1), g2=new Numero(gam2), gpl=new Numero(gampl),
gmi=new Numero(gammi); beschb(xmu,g1,g2,gpl,gmi); gam1=g1.valor; gam2=g2.valor; gampl=gpl.valor; gammi=gmi.valor; //System.out.println("bessy "+gam1); ff=2.0/PI*fact*(gam1*cosh(e)+gam2*fact2*d); e=Math.exp(e); p=e/(gampl*PI); q=1.0/(e*PI*gammi); pimu2=0.5*pimu; fact3 = (Math.abs(pimu2) < EPS ? 1.0 : Math.sin(pimu2)/pimu2); r=PI*pimu2*fact3*fact3; c=1.0; d = -x2*x2; sum=ff+r*q; sum1=p; for (i=1;i<=MAXIT;i++) { ff=(i*ff+p+q)/(i*i-xmu2); c *= (d/i); p /= (i-xmu); q /= (i+xmu); del=c*(ff+r*q); sum += del; del1=c*p-i*del; sum1 += del1; if (Math.abs(del) < (1.0+Math.abs(sum))*EPS) break; } if (i > MAXIT) nrerror("bessy series failed to converge"); rymu = -sum; ry1 = -sum1*xi2; rymup=xmu*xi*rymu-ry1; rjmu=w/(rymup-f*rymu); } else { a=0.25-xmu2; p = -0.5*xi; q=1.0; br=2.0*x; bi=2.0; fact=a*xi/(p*p+q*q); cr=br+q*fact; ci=bi+p*fact; den=br*br+bi*bi; dr=br/den; di = -bi/den; dlr=cr*dr-ci*di; dli=cr*di+ci*dr; temp=p*dlr-q*dli; q=p*dli+q*dlr; p=temp; for (i=2;i<=MAXIT;i++) { a += 2*(i-1); bi += 2.0; dr=a*dr+br; di=a*di+bi; if (Math.abs(dr)+Math.abs(di) < FPMIN) dr=FPMIN; fact=a/(cr*cr+ci*ci); cr=br+cr*fact; ci=bi-ci*fact; if (Math.abs(cr)+Math.abs(ci) < FPMIN) cr=FPMIN; den=dr*dr+di*di; dr /= den; di /= -den; dlr=cr*dr-ci*di; dli=cr*di+ci*dr; temp=p*dlr-q*dli; q=p*dli+q*dlr; p=temp; if (Math.abs(dlr-1.0)+Math.abs(dli) < EPS) break; } if (i > MAXIT) nrerror("cf2 failed in bessjy"); gam=(p-f)/q; rjmu=Math.sqrt(w/((p-f)*gam+q)); rjmu=SIGN(rjmu,rjl); rymu=rjmu*gam; rymup=rymu*(p+q/gam); ry1=xmu*xi*rymu-rymup; } fact=rjmu/rjl; rj.valor=rjl1*fact; rjp.valor=rjp1*fact; for (i=1;i<=nl;i++) { rytemp=(xmu+i)*xi2*ry1-rymu; rymu=ry1; ry1=rytemp; } ry.valor=rymu; ryp.valor=xnu*xi*rymu-ry1; } public static void bessik(double x, double xnu, Numero ri, Numero rk, Numero rip,
Numero rkp){ int i,l,nl; double a,a1,b,c,d,del,del1,delh,dels,e,f,fact,fact2,ff,gam1=0.0,gam2=0.0, gammi=0.0,gampl=0.0,h,p,pimu,q,q1,q2,qnew,ril,ril1,rimu,rip1,ripl, ritemp,rk1,rkmu,rkmup,rktemp,s,sum,sum1,x2,xi,xi2,xmu,xmu2; if (x <= 0.0 || xnu < 0.0) nrerror("bad arguments in bessik"); nl=(int)(xnu+0.5); xmu=xnu-nl; xmu2=xmu*xmu; xi=1.0/x; xi2=2.0*xi; h=xnu*xi; if (h < FPMIN) h=FPMIN; b=xi2*xnu; d=0.0; c=h; for (i=1;i<=MAXIT;i++) { b += xi2; d=1.0/(b+d); c=b+1.0/c; del=c*d; h=del*h; if (Math.abs(del-1.0) < EPS) break; } if (i > MAXIT) nrerror("x too large in bessik; try asymptotic expansion"); ril=FPMIN; ripl=h*ril; ril1=ril; rip1=ripl; fact=xnu*xi; for (l=nl;l>=1;l--) { ritemp=fact*ril+ripl; fact -= xi; ripl=fact*ritemp+ril; ril=ritemp; } f=ripl/ril; if (x < XMIN) { x2=0.5*x; pimu=PI*xmu; fact = (Math.abs(pimu) < EPS ? 1.0 : pimu/Math.sin(pimu)); d = -Math.log(x2); e=xmu*d; fact2 = (Math.abs(e) < EPS ? 1.0 : sinh(e)/e); Numero g1=new Numero(gam1), g2=new Numero(gam2), gpl=new Numero(gampl),
gmi=new Numero(gammi); beschb(xmu,g1,g2,gpl,gmi); gam1=g1.valor; gam2=g2.valor; gampl=gpl.valor; gammi=gmi.valor; ff=fact*(gam1*cosh(e)+gam2*fact2*d); sum=ff; e=Math.exp(e); p=0.5*e/gampl; q=0.5/(e*gammi); c=1.0; d=x2*x2; sum1=p; for (i=1;i<=MAXIT;i++) { ff=(i*ff+p+q)/(i*i-xmu2); c *= (d/i); p /= (i-xmu); q /= (i+xmu); del=c*ff; sum += del; del1=c*(p-i*ff); sum1 += del1; if (Math.abs(del) < Math.abs(sum)*EPS) break; } if (i > MAXIT) nrerror("bessk series failed to converge"); rkmu=sum; rk1=sum1*xi2; } else { b=2.0*(1.0+x); d=1.0/b; h=delh=d; q1=0.0; q2=1.0; a1=0.25-xmu2; q=c=a1; a = -a1; s=1.0+q*delh; for (i=2;i<=MAXIT;i++) { a -= 2*(i-1); c = -a*c/i; qnew=(q1-b*q2)/a; q1=q2; q2=qnew; q += c*qnew; b += 2.0; d=1.0/(b+a*d); delh=(b*d-1.0)*delh; h += delh; dels=q*delh; s += dels; if (Math.abs(dels/s) < EPS) break; } if (i > MAXIT) nrerror("bessik: failure to converge in cf2"); h=a1*h; rkmu=Math.sqrt(PI/(2.0*x))*Math.exp(-x)/s; rk1=rkmu*(xmu+x+0.5-h)*xi; } rkmup=xmu*xi*rkmu-rk1; rimu=xi/(f*rkmu-rkmup); ri.valor=(rimu*ril1)/ril; rip.valor=(rimu*rip1)/ril; for (i=1;i<=nl;i++) { rktemp=(xmu+i)*xi2*rk1+rkmu; rkmu=rk1; rk1=rktemp; } rk.valor=rkmu; rkp.valor=xnu*xi*rkmu-rk1; } static double airy(double x){ final double ONOVRT=0.57735027; double absx,rootx,z; double ai; absx=Math.abs(x); rootx=Math.sqrt(absx); z=2*absx*rootx/3; if (x > 0.0) { Numero ri=new Numero(0.0), rk=new Numero(0.0), rip=new Numero(0.0),
rkp=new Numero(0.0); bessik(z,1.0/3,ri,rk,rip,rkp); ai=rootx*ONOVRT*rk.valor/PI; } else if (x < 0.0) { Numero rj=new Numero(0.0), ry=new Numero(0.0), rjp=new Numero(0.0),
ryp=new Numero(0.0); bessjy(z,1.0/3,rj,ry,rjp,ryp); ai=0.5*rootx*(rj.valor-ONOVRT*ry.valor); } else { ai=0.35502805; } return ai; } static private double sinh(double x){ double y=(Math.exp(x)-Math.exp(-x))/2; return y; } static private double cosh(double x){ double y=(Math.exp(x)+Math.exp(-x))/2; return y; } static private int IMAX(int a, int b){ if(a>b) return a; return b; } static private double SIGN(double a, double b){ int sgn=1; if(b<0) sgn=-1; return (Math.abs(a)*sgn); } static void nrerror(String texto){ System.out.println(texto); } }

Llamada a la función estática airy de la clase Airy

   for(double ang =-5*Math.PI/180; ang<Math.PI/180; ang+=Math.PI/18000){
        y=funcion(ang);
//...
   }
 double funcion(double ang){ //intensidad
     double temp=4*Math.PI*a*1.0e6/(3*lonOnda);
     double z=ang*Math.sqrt(n*n-1)*Math.pow(temp, 2.0/3)/Math.pow((4-n*n), 1.0/6);
     double ai=Airy.airy(z);
    return (ai*ai);
 }
   double iRefraccion(int lonOnda){
        double z=1.4391-4.4104e-4*lonOnda+6.3095e-7*lonOnda*lonOnda-
3.1941e-10*lonOnda*lonOnda*lonOnda; return z; }

Referencia

Press W. H., Teukolsky S. A., Vetterling W. T., Flannery B. P. Numerical Recipes in C, Second edition. Special functions. Sección 6.7. Bessel functions of fractional order. Airy fucntions. Cambridge University Press. Código en C adaptado por el autor al lenguaje Java

Anterior