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

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